text:
function set_boundary(Qbase, cellxmax, cellymax, vecAx, vecAy, bdcon, Rd, g, nval)
    """bdcon[i][j]
        i : 境界番号(x-,x+ ,y-,y+)
        j=1-6 : "bd1_con":"2",
                "bd1_rho":"1.0",
                "bd1_u"  :"300.0",
                "bd1_v"  :"0.0",
                "bd1_p"  :"1.0",
                "bd1_N2" :"1.0",
                "bd1_N"  :"0.0",
                "bd1_T"  :"300.0",
    """

    # bd1 = x-
    if Int(bdcon[1][1]) == 0
        for j in 1:cellymax
            for l in 1:nval
                Qbase[1,j,l] = bdcon[1][l+1]
            end
        end
    elseif Int(bdcon[1][1]) == 1
        for j in 1:cellymax
            for l in 1:nval
                Qbase[1,j,l] = Qbase[2,j,l]
            end
        end
    elseif Int(bdcon[1][1]) == 2
        for j in 1:cellymax
            for l in 1:nval
                Qbase[1,j,l] = Qbase[2,j,l]
            end

            xvec = vecAx[2,j,1]
            yvec = vecAx[2,j,2]
            u = Qbase[2,j,2]
            v = Qbase[2,j,3]

            Qbase[1,j,2] = ((-xvec^2+yvec^2)*u-2*xvec*yvec*v)/(xvec^2+yvec^2)
            Qbase[1,j,3] = (-2*xvec*yvec*u+(xvec^2-yvec^2)*v)/(xvec^2+yvec^2)
        end
    elseif Int(bdcon[1][1]) == 3
        for j in 1:cellymax
            for l in 1:nval
                Qbase[1,j,l] = Qbase[2,j,l]
            end

            u = Qbase[2,j,2]
            v = Qbase[2,j,3]

            Qbase[1,j,2] = -u
            Qbase[1,j,3] = -v
        end
    elseif Int(bdcon[1][1]) == 4
        for j in 1:cellymax
            for l in 1:nval
                Qbase[1,j,l] = Qbase[cellxmax-1,j,l]
            end
        end
    end

    # bd2 = x+
    if Int(bdcon[2][1]) == 0
        for j in 1:cellymax
            for l in 1:nval
                Qbase[cellxmax,j,l] = bdcon[2][l+1]
            end
        end
    elseif Int(bdcon[2][1]) == 1
        for j in 1:cellymax
            for l in 1:nval
                Qbase[cellxmax,j,l] = Qbase[cellxmax-1,j,l]
            end
        end
    elseif Int(bdcon[2][1]) == 2
        for j in 1:cellymax
            for l in 1:nval
                Qbase[cellxmax,j,l] = Qbase[cellxmax-1,j,l]
            end

            xvec = vecAx[cellymax,j,1]
            yvec = vecAx[cellymax,j,2]
            u = Qbase[cellxmax-1,j,2]
            v = Qbase[cellxmax-1,j,3]

            Qbase[cellxmax,j,2] = ((-xvec^2+yvec^2)*u-2*xvec*yvec*v)/(xvec^2+yvec^2)
            Qbase[cellxmax,j,3] = (-2*xvec*yvec*u+(xvec^2-yvec^2)*v)/(xvec^2+yvec^2)
        end
    elseif Int(bdcon[2][1]) == 3
        for j in 1:cellymax
            for l in 1:nval
                Qbase[cellxmax,j,l] = Qbase[cellxmax-1,j,l]
            end
            u = Qbase[cellxmax-1,j,2]
            v = Qbase[cellxmax-1,j,3]

            Qbase[cellxmax,j,2] = -u
            Qbase[cellxmax,j,3] = -v
        end
    elseif Int(bdcon[2][1]) == 4
        for j in 1:cellymax
            for l in 1:nval
                Qbase[cellxmax,j,l] = Qbase[2,j,l]
            end
        end
    end

    # bd3 = y-
    if Int(bdcon[3][1]) == 0
        for i in 1:cellxmax
            for l in 1:nval
                Qbase[i,1,l] = bdcon[3][l+1]
            end
        end
    elseif Int(bdcon[3][1]) == 1
        for i in 1:cellxmax
            for l in 1:nval
                Qbase[i,1,l] = Qbase[i,2,l]
            end
        end
    elseif Int(bdcon[3][1]) == 2
        for i in 1:cellxmax
            for l in 1:nval
                Qbase[i,1,l] = Qbase[i,2,l]
            end

            xvec = vecAy[i,2,1]
            yvec = vecAy[i,2,2]
            u = Qbase[i,2,2]
            v = Qbase[i,2,3]

            Qbase[i,1,2] = ((-xvec^2+yvec^2)*u-2*xvec*yvec*v)/(xvec^2+yvec^2)
            Qbase[i,1,3] = (-2*xvec*yvec*u+(xvec^2-yvec^2)*v)/(xvec^2+yvec^2)
        end
    elseif Int(bdcon[3][1]) == 3
        for i in 1:cellxmax
            for l in 1:nval
                Qbase[i,1,l] = Qbase[i,2,l]
            end

            u = Qbase[i,2,2]
            v = Qbase[i,2,3]

            Qbase[i,1,2] = -u
            Qbase[i,1,3] = -v
        end
    elseif Int(bdcon[3][1]) == 4
        for i in 1:cellxmax
            for l in 1:nval
                Qbase[i,1,l] = Qbase[i,cellymax-1,l]
            end
        end
    elseif Int(bdcon[3][1]) == 7
        for i in 1:cellxmax
            for l in 1:nval
                Qbase[i,1,l] = Qbase[i,2,l]
            end

            u = Qbase[i,2,2]
            v = Qbase[i,2,3]

            Qbase[i,1,2] = -u
            Qbase[i,1,3] = -v

            Tw = bdcon[3][nval+2]
            p  = Qbase[i,1,4]
            Qbase[i,1,1] = p/(Tw*Rd)
        end
    end

    # bd4 = y+
    if Int(bdcon[4][1]) == 0
        for i in 1:cellxmax
            for l in 1:nval
                Qbase[i,cellymax,l] = bdcon[4][l+1]
            end
        end
    elseif Int(bdcon[4][1]) == 1
        for i in 1:cellxmax
            for l in 1:nval
                Qbase[i,cellymax,l] = Qbase[i,cellymax-1,l]
            end
        end
    elseif Int(bdcon[4][1]) == 2
        for i in 1:cellxmax
            for l in 1:nval
                Qbase[i,cellymax,l] = Qbase[i,cellymax-1,l]
            end

            xvec = vecAy[i,cellymax,1]
            yvec = vecAy[i,cellymax,2]
            u = Qbase[i,cellymax-1,2]
            v = Qbase[i,cellymax-1,3]

            Qbase[i,cellymax,2] = ((-xvec^2+yvec^2)*u-2*xvec*yvec*v)/(xvec^2+yvec^2)
            Qbase[i,cellymax,3] = (-2*xvec*yvec*u+(xvec^2-yvec^2)*v)/(xvec^2+yvec^2)
        end
    elseif Int(bdcon[4][1]) == 3
        for i in 1:cellxmax
            for l in 1:nval
                Qbase[i,cellymax,l] = Qbase[i,cellymax-1,l]
            end

            u = Qbase[i,cellymax-1,2]
            v = Qbase[i,cellymax-1,3]

            Qbase[i,cellymax,2] = -u
            Qbase[i,cellymax,3] = -v
        end
    elseif Int(bdcon[4][1]) == 4
        for i in 1:cellxmax
            for l in 1:nval
                Qbase[i,cellymax,l] = Qbase[i,2,l]
            end
        end
    elseif Int(bdcon[4][1]) == 5
        temp1 = Int64(round(cellxmax/4)+1) # 1/4流出
        temp2 = Int64((temp1-1)*3+1)       # 1/4~3/4流入
        for i in 1:temp1
            for l in 1:nval
                Qbase[i,cellymax,l] = Qbase[i,cellymax-1,l]
            end
        end
        for i in temp1+1:temp2
            for l in 1:nval
                Qbase[i,cellymax,l] = bdcon[4][l+1]
            end
        end
        for i in temp2+1:cellxmax
            for l in 1:nval
                Qbase[i,cellymax,l] = Qbase[i,cellymax-1,l]
            end
        end
    elseif Int(bdcon[4][1]) == 6
        temp1 = Int64(round(cellxmax/4)+1) # 1/4流出
        temp2 = Int64((temp1-1)*3+1)       # 1/4~3/4流入
        for i in 1:temp1
            for l in 1:nval
                Qbase[i,cellymax,l] = Qbase[i,cellymax-1,l]
            end
        end
        for i in temp1+1:temp2
            for l in 1:nval
                Qbase[i,cellymax,l] = bdcon[4][l+1]
            end
            T = bdcon[4][nval+2]
            p = (Qbase[i,cellymax,1]*Rd) * T
            Qbase[i,cellymax,4] = p
        end
        for i in temp2+1:cellxmax
            for l in 1:nval
                Qbase[i,cellymax,l] = Qbase[i,cellymax-1,l]
            end
        end
    end

    return Qbase
end

function check_bd(bdcon)
    for l in 1:4
        if Int(bdcon[l][1]) == 0
        elseif Int(bdcon[l][1]) == 1
        elseif Int(bdcon[l][1]) == 2
        elseif Int(bdcon[l][1]) == 3
        elseif Int(bdcon[l][1]) == 4
        elseif Int(bdcon[l][1]) == 5
        elseif Int(bdcon[l][1]) == 6
        elseif Int(bdcon[l][1]) == 7
        else
            println("\n check boundary condition ! \n")
            throw(UndefVarError(:x))
        end
    end
end


text:
using Gtk

win = GtkWindow("My First Gtk.jl Program", 400, 200)

b = GtkButton("Click Me")
push!(win,b)

function on_button_clicked(w)
    println("The button has been clicked")
end

signal_connect(on_button_clicked, b, "clicked")

showall(win)

title = get_gtk_property(win, :title, String)

ext: jl


text:
function Table_SINDY_PI(s, epsilon, npoly, ntrig, w)

    nvar = length(s)
    table = []
    #push!(table, "")

    push!(table, "1")

   if npoly >=1
    for i in 1:nvar
        push!(table, s[i])
    end
end

    if npoly >= 2
        for i in 1:nvar
            for j in i:nvar
                push!(table, string(s[i], s[j]))
            end
        end
    end

    if npoly >= 3
        for i in 1:nvar
            for j in i:nvar
                for k in j:nvar
                  push!(table, string(s[i], s[j], s[k]))
            end
            end
        end
    end

    if npoly >= 4
        for i in 1:nvar
            for j in i:nvar
                for k in j:nvar
                    for l in k:nvar
                      push!(table, string(s[i], s[j], s[k], s[l]))

            end
            end
        end
    end
end

if npoly >= 5
    for i in 1:nvar
        for j in i:nvar
            for k in j:nvar
                for l in k:nvar
                    for m in l:nvar
                  push!(table, string(s[i], s[j], s[k], s[l], s[m]))

        end
        end
    end
end
end
end

if npoly >= 6
    for i in 1:nvar
        for j in i:nvar
            for k in j:nvar
                for l in k:nvar
                    for m in l:nvar
                        for n in m:nvar
                  push!(table, string(s[i], s[j], s[k], s[l], s[m], s[n]))

        end
        end
    end
end
end
end
end

if ntrig >=1
   for i in 1:nvar
       push!(table, string("sin", "", s[i]))
   end

   for i in 1:nvar
       push!(table, string("cos", "", s[i]))
   end
end

if ntrig >=2
   for i in 1:nvar
       push!(table, string("cos^{2}", "", s[i]))
   end
end

if ntrig >=3
    for i in 1:nvar
        for j in 1:nvar
            push!(table, string(s[i], "*", s[i], "sin", "", s[j], "cos", "", s[j]))
        end
    end
end

if ntrig >=4
    for i in 1:nvar
        for j in 1:nvar
            push!(table, string(s[i], "*", s[i], "sin", "", s[j]))
        end
    end

    for i in 1:nvar
        for j in 1:nvar
            push!(table, string(s[i], "*", s[i], "cos", "", s[j]))
        end
    end

    for i in 1:nvar
        for j in 1:nvar
            push!(table, string("sin", "", s[i], "cos", "", s[j]))
        end
    end
end

tablecopy = deepcopy(table)

for i in 1:length(table)-1
    push!(table, string("dot", "(", s[w], ")" , "", tablecopy[i+1]))
end

table = hcat(table, epsilon)

    return table
end


text:
module NewtonsMethod

greet() = print("Hello World!")

using LinearAlgebra
function newtonroot(f, f_prime, x_0; tolerance= 1E-7, maxiter= 1000)
    x_old = x_0
    error = Inf
    iter = 1
    while iter <= maxiter && error > tolerance
        x_new = x_old - (f(x_old)/f_prime(x_old))
        error = norm(x_new - x_old)
        x_old = x_new
        iter  = iter+1
    end
    if error > tolerance
        return nothing
    else
    return (x_old)
    end
end

using ForwardDiff

D(f) = x -> ForwardDiff.derivative(f, x)

function newtonroot(f, x_0; tolerance= 1E-7, maxiter= 1000)
    x_old = x_0
    error = Inf
    iter = 1
    f_prime =  D(f)
    while iter <= maxiter && error > tolerance
        x_new = x_old - (f(x_old)/f_prime(x_old))
        error = norm(x_new - x_old)
        x_old = x_new
        iter  = iter+1
    end
    if error > tolerance
    return nothing
    else
    return (x_old)
    end
end

export newtonroot
end # module


text:
# generate fake XSV data

function genxsv(;num_rows  = 100,
                 num_cols  = 10,
                 quotechar = '"',
                 delimiter = ',',
                 max_field = 20,
                 new_line  = "\n")
	field_data = Array(ASCIIString, 0)
    row_data   = Array(Vector{ASCIIString}, 0)
    rows       = Array(Vector{ASCIIString}, 0)

    for i=1:num_rows
        fields     = Array(ASCIIString, 0)
        field_data = Array(ASCIIString, 0)
        for i=1:num_cols
            field = ASCIIString(rand(0x00:0x7f, rand(1:max_field)))
            push!(field_data, field)
            if in(quotechar, field)
            	field = join(split(field, quotechar), string(quotechar)*string(quotechar))
            end
            if in(quotechar, field) || in(delimiter, field) || contains(new_line, field) || in('\n', field) || in('\r', field)
                field = string(quotechar) * field * string(quotechar)
            end
            push!(fields, field)
        end
        push!(row_data, field_data)
        push!(rows, fields)
    end
    join([join(row, delimiter) for row in rows], new_line), row_data
end


text:
# Physical constants

const cp = 1005          # Specific heat capacity of dry air (J/K/kg)
const epsm = 0.622       # Ratio of molecular weights of water and dry air    BEFORE eps
const e0 = 611.213       # Saturation vapour pressure at Tm (Pa)
const g = 9.81           # Acceleration due to gravity (m/s^2)
const hcap_ice = 2100.   # Specific heat capacity of ice (J/K/kg)
const hcap_wat = 4180.   # Specific heat capacity of water (J/K/kg)
const hcon_air = 0.025   # Thermal conductivity of air (W/m/K)
const hcon_clay = 1.16   # Thermal conductivity of clay (W/m/K)
const hcon_ice = 2.24    # Thermal conducivity of ice (W/m/K)
const hcon_sand = 1.57   # Thermal conductivity of sand (W/m/K)
const hcon_wat = 0.56    # Thermal conductivity of water (W/m/K)
const Lc = 2.501e6       # Latent heat of condensation (J/kg)
const Lf = 0.334e6       # Latent heat of fusion (J/kg)
const Ls = Lc + Lf       # Latent heat of sublimation (J/kg)
const R = 8.3145         # Universal gas constant (J/K/mol)
const Rgas = 287         # Gas constant for dry air (J/K/kg)
const Rwat = 462         # Gas constant for water vapour (J/K/kg)
const rho_ice = 917.     # Density of ice (kg/m^3)
const rho_wat = 1000.    # Density of water (kg/m^3)
const sb = 5.67e-8       # Stefan-Boltzmann constant (W/m^2/K^4)
const Tm = 273.15        # Melting point (K)
const vkman = 0.4        # Von Karman constant


text:
# Error type and decoding functionality

export CuError

"""
    CuError(code)
    CuError(code, meta)

Create a CUDA error object with error code `code`. The optional `meta` parameter indicates
whether extra information, such as error logs, is known.
"""
struct CuError <: Exception
    code::CUresult
    meta::Any

    CuError(code, meta=nothing) = new(code, meta)
end

Base.convert(::Type{CUresult}, err::CuError) = err.code

Base.:(==)(x::CuError,y::CuError) = x.code == y.code

"""
    name(err::CuError)

Gets the string representation of an error code.

```jldoctest
julia> err = CuError(CUDA.cudaError_enum(1))
CuError(CUDA_ERROR_INVALID_VALUE)

julia> name(err)
"ERROR_INVALID_VALUE"
```
"""
function name(err::CuError)
    str_ref = Ref{Cstring}()
    cuGetErrorName(err, str_ref)
    unsafe_string(str_ref[])[6:end]
end

"""
    description(err::CuError)

Gets the string description of an error code.
"""
function description(err::CuError)
    if err.code == -1%UInt32
        "Cannot use the CUDA stub libraries."
    else
        str_ref = Ref{Cstring}()
        cuGetErrorString(err, str_ref)
        unsafe_string(str_ref[])
    end
end

function Base.showerror(io::IO, err::CuError)
    try
        print(io, "CUDA error: $(description(err)) (code $(reinterpret(Int32, err.code)), $(name(err)))")
    catch
        # we might throw before the library is initialized
        print(io, "CUDA error (code $(reinterpret(Int32, err.code)), $(err.code))")
    end

    if err.meta != nothing
        print(io, "\n")
        print(io, err.meta)
    end
end

Base.show(io::IO, ::MIME"text/plain", err::CuError) = print(io, "CuError($(err.code))")

@enum_without_prefix cudaError_enum CUDA_

## API call wrapper

# outlined functionality to avoid GC frame allocation
@noinline function initialize_api()
    prepare_cuda_state()
    return
end
@noinline function throw_api_error(res)
    if res == ERROR_OUT_OF_MEMORY
        throw(OutOfGPUMemoryError())
    else
        throw(CuError(res))
    end
end

macro check(ex)
    quote
        res = $(esc(ex))
        if res != SUCCESS
            throw_api_error(res)
        end

        nothing
    end
end


text:
using DataDeps

register(DataDep("English WikiCorpus v1.0",
    """
    Website: http://www.cs.upc.edu/~nlp/wikicorpus/
    Author: Samual Reese, et al.
    The Wikicorpus contains large portions of the Wikipedia, based on a 2006 dump.
    This is the English portion which has around 600 million words.

    Please cite the following publication if you use the corpora:
    Samuel Reese, Gemma Boleda, Montse Cuadros, Lluís Padró, German Rigau. Wikicorpus: A Word-Sense Disambiguated Multilingual Wikipedia Corpus. In Proceedings of 7th Language Resources and Evaluation Conference (LREC'10), La Valleta, Malta. May, 2010.
    """,
    "http://www.cs.upc.edu/~nlp/wikicorpus/raw.en.tgz",
    "cab102a59f7e5bfb4ac832c15210bd22840bbed5508899f2dee04b857830e858";
    post_fetch_method = unpack
))


text:
function dayNumber(inDate::Date)
	return Dates.date2epochdays(inDate)::Integer;
end

excelcostant=693959;

"""
From Excel Number Format to Date

		Date=fromExcelNumberToDate(ExcelNumber)

Where:\n
		ExcelNumber = Integer representing a date in the excel format.

		Date      = Date representing the input in the Julia object format.

# Example
```julia-repl
julia> fromExcelNumberToDate(45000)
2023-03-15
```
"""
function fromExcelNumberToDate(num::Integer)::Date
	return Dates.Dates.epochdays2date(excelcostant+num);
end

"""
Actual Number of days between two dates

		ndays=daysact(Date1,Date2)

Where:\n
		Date1 = Start date.
		Date2 = End date.

		ndays      = Actual Number of days between Start Date and End Date.

# Example
```julia-repl
julia> daysact(Date(1996,10,12),Date(1998,1,10))
455
```
"""
function daysact(Date1::Date,Date2::Date)::Integer
	D1= dayNumber(Date1);
	D2= dayNumber(Date2);
	dayCount=D2-D1;
	return dayCount;
end

function isLastOfFebruary(inDate::Date)::Bool
	return (Dates.isleapyear(Dates.year(inDate))&&(Dates.day(inDate)==29)&&(Dates.month(inDate)==2))||((!Dates.isleapyear(Dates.year(inDate)))&&(Dates.day(inDate)==28)&&(Dates.month(inDate)==2));
end

currMaxImplemented=7;

"""
Fraction of year between two Dates according the following convention

		yfr=yearfrac(Date1,Date2,basis)

Where:\n
		Date1 = Start date.
		Date2 = End date.
		basis = Integer representing the following conventions:
				- 0 = (ACT/ACT)
				- 1 = (30/360 SIA)
				- 2 = (ACT/360)
				- 3 = (ACT/365)
				- 4 = (30/360 PSA)
				- 5 = (30/360 ISDA)
				- 6 = (30E/360)
				- 7 = (ACT/365 JPN)

		yfr      = fraction of year between start and end date according to basis.

# Example
```julia-repl
julia> yearfrac(Date(1996,10,12),Date(1998,1,10),1)
1.2444444444444445
```
"""
function yearfrac(startDate::Date,endDate::Date,convention::Integer)::Real
if (convention<0)
	error("Negative basis are not defined, check the help")
end
if(convention>currMaxImplemented)
	error("Convention not implemented yet")
end
if(startDate>endDate)
	return -yearfrac(endDate,startDate,convention);
elseif(startDate==endDate)
	return 0.0;
end
yearFrac=0.0;
y1=Dates.year(startDate);
m1=Dates.month(startDate);
y2=Dates.year(endDate);
m2=Dates.month(endDate);
d1=Dates.day(startDate);
d2=Dates.day(endDate);
if (convention==0)#(ACT/ACT)
	Nday=daysact(startDate,endDate);
	EndOFYear=Date(1);
	if(isLastOfFebruary(startDate)&&(Dates.day(startDate)==29))
		EndOFYear=Date(Dates.year(startDate)+1,3,1);
	else
		EndOFYear=Date(Dates.year(startDate)+1,Dates.month(startDate),Dates.day(startDate));
	end
	yearFrac=(Nday)/ daysact(startDate,EndOFYear);
elseif(convention==1)  	#(30/360 SIA)
	if(isLastOfFebruary(startDate)&&isLastOfFebruary(endDate))
		d2=30;
	end
	if(isLastOfFebruary(startDate)||Dates.day(startDate)==31)
		d1=30;
	end
	if(d1==30&&d2==31)
		d2=30;
	end
	dy=y2-y1;
	dm=m2-m1;
	dd=d2-d1;
	yearFrac=(360.0*dy+30.0*dm+dd)/360.0;
elseif(convention==2)		#(ACT/360)
	Nday=daysact(startDate,endDate);
	yearFrac=(Nday)/360.0;
elseif(convention==3)		# (ACT/365)
	Nday=daysact(startDate,endDate);
	yearFrac=( Nday)/365.0;
elseif(convention==4)		# (30/360 PSA)
	if((Dates.day(startDate)==31)||isLastOfFebruary(startDate))
		d1=30;
	end
	if((Dates.day(startDate)==30||isLastOfFebruary(startDate))&&Dates.day(endDate)==31)
		d2=30;
	end
	dy=y2-y1;
	dm=m2-m1;
	dd=d2-d1;
	yearFrac=(360.0*dy+30.0*dm+dd)/360.0;

elseif(convention==5)		#(30/360 ISDA)
	y1=Dates.year(startDate);
	m1=Dates.month(startDate);
	y2=Dates.year(endDate);
	m2=Dates.month(endDate);
	if(Dates.day(startDate)<31)
		d1=Dates.day(startDate);
	else
		d1=30;
	end
	if ((Dates.day(endDate)==31)&&(d1>29))
		d2=30;
	else
		d2=Dates.day(endDate);
	end
	yearFrac=(360.0*((y2-y1))+30.0*((m2-m1))+(d2-d1))/360.0;

elseif(convention==6)		#(30E/360)
	if(d1==31)
		d1=30;
	end
	if(d2==31)
		d2=30;
	end
	dy=y2-y1;
	dm=m2-m1;
	dd=d2-d1;
	yearFrac=(360.0*dy+30.0*dm+dd)/360.0;
elseif(convention==7)		#(ACT/365 JPN)
	daydistance = [0;31;59;90;120;151;181;212;243;273;304;334];
	dy=y2-y1;
	dd=d2-d1;
	dayCount=(365.0*dy+daydistance[m2]-daydistance[m1]+dd)
	yearFrac = dayCount/365.0;
end

return yearFrac;

end


text:
using DashBootstrapComponents

input_groups = html_div([
    dbc_inputgroup(
        [dbc_inputgrouptext("@"), dbc_input(placeholder = "Username")],
        className = "mb-3",
    ),
    dbc_inputgroup(
        [
            dbc_input(placeholder = "Recipient's username"),
            dbc_inputgrouptext("@example.com"),
        ],
        className = "mb-3",
    ),
    dbc_inputgroup(
        [
            dbc_inputgrouptext("\$"),
            dbc_input(placeholder = "Amount", type = "number"),
            dbc_inputgrouptext(".00"),
        ],
        className = "mb-3",
    ),
    dbc_inputgroup(
        [
            dbc_inputgrouptext("Total:"),
            dbc_inputgrouptext("\$"),
            dbc_input(placeholder = "Amount", type = "number"),
            dbc_inputgrouptext(".00"),
            dbc_inputgrouptext("only"),
        ],
        className = "mb-3",
    ),
    dbc_inputgroup(
        [dbc_inputgrouptext("With textarea"), dbc_textarea()],
        className = "mb-3",
    ),
    dbc_inputgroup([
        dbc_select(
            options = [
                Dict("label" => "Option 1", "value" => 1),
                Dict("label" => "Option 2", "value" => 2),
            ],
        ),
        dbc_inputgrouptext("With select"),
    ]),
]);

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
"""

# cmdline

Recursively parse the model to construct command line.

### Method
```julia
cmdline(m)
```

### Required arguments
```julia
* `m::CmdStanSampleModel`                : CmdStanSampleModel
```

### Related help
```julia
?CmdStanSampleModel                      : Create a CmdStanSampleModel
```
"""
function cmdline(m::Union{OptimizeModel, Optimize, Lbfgs, Bfgs, Newton}, id)

  #=
  `/Users/rob/.julia/dev/StanOptimize/examples/Bernoulli/tmp/bernoulli
  optimize algorithm=lbfgs init_alpha=0.001 tol_obj=1.0e-8 tol_rel_
  obj=10000.0 tol_grad=1.0e-8 tol_rel_grad=1.0e7 tol_param=1.0e-8
  history_size=5 iter=2000 save_iterations=1 random seed=-1 init=2
  id=1 data file=/Users/rob/.julia/dev/StanOptimize/examples/Bernoulli/tmp/bernoulli_data_1.R
  output file=/Users/rob/.julia/dev/StanOptimize/examples/Bernoulli/tmp/bernoulli_chain_1.csv
  refresh=100`
  =#

  cmd = ``
  if isa(m, OptimizeModel)
    # Handle the model name field for unix and windows
    cmd = `$(m.exec_path)`

    # Sample() specific portion of the model
    cmd = `$cmd $(cmdline(getfield(m, :method), id))`

    # Common to all models
    cmd = `$cmd random seed=$(getfield(m, :seed).seed)`

    # Init file required?
    if length(m.init_file) > 0 && isfile(m.init_file[id])
      cmd = `$cmd init=$(m.init_file[id])`
    else
      cmd = `$cmd init=$(m.init.bound)`
    end

    # Data file required?
    if length(m.data_file) > 0 && isfile(m.data_file[id])
      cmd = `$cmd id=$(id) data file=$(m.data_file[id])`
    end

    # Output options
    cmd = `$cmd output`
    if length(getfield(m, :output).file) > 0
      cmd = `$cmd file=$(string(getfield(m, :output).file))`
    end
    if length(m.diagnostic_file) > 0
      cmd = `$cmd diagnostic_file=$(string(getfield(m, :output).diagnostic_file))`
    end
    cmd = `$cmd refresh=$(string(getfield(m, :output).refresh))`

  else

    # The 'recursive' part
    if isa(m, OptimizeAlgorithm)
      cmd = `$cmd algorithm=$(split(lowercase(string(typeof(m))), '.')[end])`
    else
      cmd = `$cmd $(split(lowercase(string(typeof(m))), '.')[end])`
    end
    for name in fieldnames(typeof(m))
      if  isa(getfield(m, name), String) || isa(getfield(m, name), Tuple)
        cmd = `$cmd $(name)=$(getfield(m, name))`
      elseif length(fieldnames(typeof(getfield(m, name)))) == 0
        if isa(getfield(m, name), Bool)
          cmd = `$cmd $(name)=$(getfield(m, name) ? 1 : 0)`
        else
          if name == :metric || isa(getfield(m, name), DataType)
            cmd = `$cmd $(name)=$(split(lowercase(string(typeof(getfield(m, name)))), '.')[end])`
          else
            cmd = `$cmd $(name)=$(getfield(m, name))`
          end
        end
      else
        cmd = `$cmd $(cmdline(getfield(m, name), id))`
      end
    end
  end

  cmd

end


text:
SubWorker{T,A,S,PT} = RemoteChannel{Channel{Vector{SubProblem{T,A,S,PT}}}}
ScenarioProblemChannel{S} = RemoteChannel{Channel{ScenarioProblems{S}}}

Work = RemoteChannel{Channel{Int}}
Progress{T <: AbstractFloat} = Tuple{Int,Int,T}
ProgressQueue{T <: AbstractFloat} = RemoteChannel{Channel{Progress{T}}}

function initialize_subproblems!(ph::AbstractProgressiveHedgingSolver, subsolver::QPSolver, penaltyterm::PenaltyTerm, subworkers::Vector{SubWorker{T,A,S,PT}}) where {T <: AbstractFloat, A <: AbstractVector, S <: LQSolver, PT <: PenaltyTerm}
    # Create subproblems on worker processes
    m = ph.stochasticprogram
    @sync begin
        for w in workers()
            subworkers[w-1] = RemoteChannel(() -> Channel{Vector{SubProblem{T,A,S,PT}}}(1), w)
            @async load_worker!(scenarioproblems(m), m, w, subworkers[w-1], subsolver, penaltyterm)
        end
        # Prepare memory
        log_val = ph.parameters.log
        ph.parameters.log = false
        log!(ph)
        ph.parameters.log = log_val
    end
    update_iterate!(ph)
    return ph
end

function update_dual_gap!(ph::AbstractProgressiveHedgingSolver, subworkers::Vector{SubWorker{T,A,S,PT}}) where {T <: AbstractFloat, A <: AbstractVector, S <: LQSolver, PT <: PenaltyTerm}
    # Update δ₂
    partial_δs = Vector{Float64}(undef, nworkers())
    @sync begin
        for (i,w) in enumerate(workers())
            @async partial_δs[i] = remotecall_fetch((sw,ξ)->begin
                subproblems = fetch(sw)
                if length(subproblems) > 0
                    return sum([s.π*norm(s.x-ξ,2)^2 for s in subproblems])
                else
                    return zero(T)
                end
            end,
            w,
            subworkers[w-1],
            ph.ξ)
        end
    end
    ph.data.δ₂ = sum(partial_δs)
    return nothing
end

function calculate_objective_value(ph::AbstractProgressiveHedgingSolver, subworkers::Vector{<:SubWorker{T}}) where T <: AbstractFloat
    partial_objectives = Vector{Float64}(undef, nworkers())
    @sync begin
        for (i,w) in enumerate(workers())
            @async partial_objectives[i] = remotecall_fetch((sw)->begin
                subproblems = fetch(sw)
                if length(subproblems) > 0
                    return sum([get_objective_value(s) for s in subproblems])
                else
                    return zero(T)
                end
            end,
            w,
            subworkers[w-1])
        end
    end
    return sum(partial_objectives)
end

function fill_submodels!(ph::AbstractProgressiveHedgingSolver, scenarioproblems::ScenarioProblems, subworkers::Vector{<:SubWorker}, nrows::Integer, ncols::Integer)
    j = 0
    @sync begin
        for w in workers()
            n = remotecall_fetch((sw)->length(fetch(sw)), w, subworkers[w-1])
            for i = 1:n
                k = i+j
                @async fill_submodel!(scenarioproblems.problems[k],remotecall_fetch((sw,i,x)->begin
                    sp = fetch(sw)[i]
                    get_solution(sp)
                end,
                w,
                subworkers[w-1],
                i,
                ph.ξ)..., nrows, ncols)
            end
            j += n
        end
    end
end

function fill_first_stage!(ph::AbstractProgressiveHedgingSolver, stochasticprogram::StochasticProgram, subworkers::Vector{<:SubWorker}, nrows::Integer, ncols::Integer)
    μ, λ = remotecall_fetch((sw,nrows,ncols) -> begin
        subproblem = fetch(sw)[1]
        μ = try
            MPB.getreducedcosts(subproblem.solver.lqmodel)[1:ncols]
        catch
            fill(NaN, ncols)
        end
        λ = try
            MPB.getconstrduals(subproblem.solver.lqmodel)[1:nrows]
        catch
            fill(NaN, nrows)
        end
        μ, λ
    end,
    2,
    subworkers[1],
    nrows,
    ncols)
    StochasticPrograms.set_first_stage_redcosts!(stochasticprogram, μ)
    StochasticPrograms.set_first_stage_duals!(stochasticprogram, λ)
    return nothing
end

function fill_submodels!(ph::AbstractProgressiveHedgingSolver, scenarioproblems::DScenarioProblems, subworkers::Vector{<:SubWorker}, nrows::Integer, ncols::Integer)
    @sync begin
        for w in workers()
            @async remotecall(fill_submodels!,
                              w,
                              subworkers[w-1],
                              ph.ξ,
                              scenarioproblems[w-1],
                              nrows,
                              ncols)
        end
    end
end

function load_worker!(scenarioproblems::ScenarioProblems,
                      sp::StochasticProgram,
                      w::Integer,
                      worker::SubWorker,
                      subsolver::QPSolver,
                      penaltyterm::PenaltyTerm)
    n = StochasticPrograms.nscenarios(scenarioproblems)
    (nscen, extra) = divrem(n, nworkers())
    prev = [nscen + (extra + 2 - p > 0) for p in 2:(w-1)]
    start = isempty(prev) ? 1 : sum(prev) + 1
    stop = min(start + nscen + (extra + 2 - w > 0) - 1, n)
    return remotecall_fetch(initialize_subworker!,
                            w,
                            worker,
                            generator(sp, :stage_1),
                            generator(sp, :stage_2),
                            stage_parameters(sp, 1),
                            stage_parameters(sp, 2),
                            scenarios(sp)[start:stop],
                            decision_length(sp),
                            subsolver,
                            penaltyterm,
                            start)
end

function load_worker!(scenarioproblems::DScenarioProblems,
                      sp::StochasticProgram,
                      w::Integer,
                      worker::SubWorker,
                      subsolver::QPSolver,
                      penaltyterm::PenaltyTerm)
    leading_scen = [scenarioproblems.scenario_distribution[p-1] for p in 2:(w-1)]
    start_id = isempty(leading_scen) ? 1 : sum(leading_scen)+1
    return remotecall_fetch(initialize_subworker!,
                            w,
                            worker,
                            generator(sp, :stage_1),
                            generator(sp, :stage_2),
                            stage_parameters(sp, 1),
                            stage_parameters(sp, 2),
                            scenarioproblems[w-1],
                            decision_length(sp),
                            subsolver,
                            penaltyterm,
                            start_id)
end

function initialize_subworker!(subworker::SubWorker{T,A,S,PT},
                               stage_one_generator::Function,
                               stage_two_generator::Function,
                               stage_one_params::Any,
                               stage_two_params::Any,
                               scenarios::Vector{<:AbstractScenario},
                               xdim::Integer,
                               subsolver::QPSolver,
                               penaltyterm::PenaltyTerm,
                               start_id::Integer) where {T <: AbstractFloat, A <: AbstractArray, S <: LQSolver, PT <: PenaltyTerm}
    subproblems = Vector{SubProblem{T,A,S,PT}}(undef, length(scenarios))
    id = start_id
    solver = get_solver(subsolver)
    for (i,scenario) = enumerate(scenarios)
        subproblems[i] = SubProblem(_WS(stage_one_generator, stage_two_generator, stage_one_params, stage_two_params, scenario, solver),
                                    id,
                                    probability(scenario),
                                    xdim,
                                    solver,
                                    penaltyterm)
        id += 1
    end
    put!(subworker, subproblems)
end

function initialize_subworker!(subworker::SubWorker{T,A,S,PT},
                               stage_one_generator::Function,
                               stage_two_generator::Function,
                               stage_one_params::Any,
                               stage_two_params::Any,
                               scenarioproblems::ScenarioProblemChannel,
                               xdim::Integer,
                               subsolver::QPSolver,
                               penaltyterm::PenaltyTerm,
                               start_id::Integer) where {T <: AbstractFloat, A <: AbstractArray, S <: LQSolver, PT <: PenaltyTerm}
    sp = fetch(scenarioproblems)
    subproblems = Vector{SubProblem{T,A,S,PT}}(undef, StochasticPrograms.nscenarios(sp))
    id = start_id
    solver = get_solver(subsolver)
    for (i,scenario) = enumerate(scenarios(sp))
        subproblems[i] = SubProblem(_WS(stage_one_generator, stage_two_generator, stage_one_params, stage_two_params, scenario, solver),
                                    id,
                                    probability(scenario),
                                    xdim,
                                    solver,
                                    penaltyterm)
        id += 1
    end
    put!(subworker, subproblems)
end

function fill_submodels!(subworker::SubWorker{T,A,S,PT},
                         x::A,
                         scenarioproblems::ScenarioProblemChannel,
                         nrows::Integer,
                         ncols::Integer) where {T <: AbstractFloat, A <: AbstractArray, S <: LQSolver, PT <: PenaltyTerm}
    sp = fetch(scenarioproblems)
    subproblems::Vector{SubProblem{T,A,S,PT}} = fetch(subworker)
    for (i, submodel) in enumerate(sp.problems)
        fill_submodel!(submodel, subproblems[i], nrows, ncols)
    end
end

function resolve_subproblems!(subworker::SubWorker{T,A,S,PT}, ξ::AbstractVector, r::AbstractFloat) where {T <: AbstractFloat, A <: AbstractArray, S <: LQSolver, PT <: PenaltyTerm}
    subproblems::Vector{SubProblem{T,A,S,PT}} = fetch(subworker)
    Qs = A(undef, length(subproblems))
    for (i,subproblem) ∈ enumerate(subproblems)
        reformulate_subproblem!(subproblem, ξ, r)
        Qs[i] = subproblem(ξ)
    end
    return sum(Qs)
end

function collect_primals(subworker::SubWorker{T,A,S,PT}, n::Integer) where {T <: AbstractFloat, A <: AbstractArray, S <: LQSolver, PT <: PenaltyTerm}
    subproblems::Vector{SubProblem{T,A,S,PT}} = fetch(subworker)
    if length(subproblems) > 0
        return sum([subproblem.π*subproblem.x for subproblem in subproblems])
    else
        return zeros(T,n)
    end
end


text:
using Distributed

@everywhere using Pkg
@everywhere Pkg.activate("./")
@everywhere using JLD2, FileIO
@everywhere using Random
@everywhere using RandomNumbers.Xorshifts
@everywhere using ProgressMeter
@everywhere using Printf
@everywhere using Boltzmann
# using Base.Threads

# Total number of confs
freq = 1_000
total_sweeps = 10^9
nconfs = round(Int, total_sweeps / freq)
Ts = 1.2:0.1:3.4
# System size
L = 8
# Generate the Ising configurations
@everywhere include(joinpath("scripts", "ising.jl"))
pmap(x -> ising_simple(freq, L, total_sweeps, x), Ts)
# Now fit the RBMs
@everywhere include(joinpath("scripts", "fit.jl"))
pmap(x -> fit_ising(x; L=L, nconfs=nconfs), Ts)

# We now analyse everything
include("analyse.jl")
compute_values(;nconfs=nconfs, L=L)

max_stars_repo_name: edwinb-ai/IsingRBM

max_issues_repo_name: edwinb-ai/IsingRBM

max_forks_repo_name: edwinb-ai/IsingRBM


text:
module Filtering

using DataFrames, Query

######################
## single condition ##
######################

# Query - stand alone
sweet_tooth_sa = person |>
    @filter(_.sweet_tooth=="yes") |>
    DataFrame

# Query - LINQ
sweet_tooth_q = @from i in person begin
    @where i.sweet_tooth=="yes"
    @select i
    @collect DataFrame
end

#########################
## multiple conditions ##
#########################

# Query - stand alone
sweet_tooth_older_sa = person |>
    @filter(_.sweet_tooth=="yes" && _.age > 20) |>
    DataFrame

# Query - LINQ
sweet_tooth_older_q = @from i in person begin
    @where i.sweet_tooth=="yes" && i.age > 20
    @select i
    @collect DataFrame
end

end # module

max_stars_repo_name: curtisalexander/saspyrilia-sandboxes

max_issues_repo_name: curtisalexander/saspyrilia-sandboxes

max_forks_repo_name: curtisalexander/saspyrilia-sandboxes


text:
function tuning_display(p)
    lines=show(p.output, p.tuner, progress=true)
    println();
    return lines
end

function convergence_display(p)
    try
        ciplot=lineplot([p.counter-(length(p.convergence_history)-1):p.counter...], p.convergence_history, title="Convergence Interval Recent History", xlabel="Iterate",ylabel="CI", color=:yellow)
        lines=nrows(ciplot.graphics)+5
        show(p.output, ciplot); println()
        return lines
    catch
        printstyled(p.output, "CONVERGENCE PLOT UNAVAILABLE. STANDBY\n", bold=true, color=:yellow); println()
        return 2
    end
end

function evidence_display(p)
    try
        log_Zis=p.e.log_Zi[2:end]
        posidx=findfirst(i->i > -Inf, log_Zis)
        log_Zis[1:posidx-1].=log_Zis[posidx]
        evplot=lineplot(log_Zis, title="Evidence History", xlabel="Iterate", color=:red, name="Ensemble logZ")
        lines=nrows(evplot.graphics)+5
        show(p.output, evplot); println()
        return lines
    catch
        printstyled(p.output, "EVIDENCE PLOT UNAVAILABLE. STANDBY\n", bold=true, color=:red); println()
        return 2
    end
end

function info_display(p)
    try
        infoplot=lineplot(p.e.Hi[2:end], title="Information History", xlabel="Iterate", color=:green, name="Ensemble H")
        lines=nrows(infoplot.graphics)+5
        show(p.output, infoplot); println()
        return lines
    catch
        printstyled(p.output, "INFORMATION PLOT UNAVAILABLE. STANDBY\n", bold=true, color=:green); println()
        return 2
    end
end

function lh_display(p)
    try
        lhplot=lineplot(p.e.log_Li[2:end], title="Contour History", xlabel="Iterate", color=:magenta, name="Ensemble logLH")
        lines=nrows(lhplot.graphics)+5
        show(p.output, lhplot); println()
        return lines
    catch
        printstyled(p.output, "CONTOUR HISTORY UNAVAILABLE. STANDBY\n", bold=true, color=:magenta); println()
        return 2
    end
end

function liwi_display(p)
    try
        liwiplot=lineplot([max(2,p.counter-(CONVERGENCE_MEMORY-1)):p.counter...],p.e.log_Liwi[max(2,end-(CONVERGENCE_MEMORY-1)):end], title="Recent iterate evidentiary weight", xlabel="Iterate", name="Ensemble log Liwi", color=:cyan)
        lines=nrows(liwiplot.graphics)+5
        show(p.output, liwiplot); println()
        return lines
    catch
        printstyled(p.output, "EVIDENTIARY HISTORY UNAVAILABLE. STANDBY\n", bold=true, color=:cyan); println()
        return 2
    end
end

function ensemble_display(p)
    return lines=show(p.output, p.e, progress=true)
end

function model_display(p)
    try
        println("Current MAP model:")
        return lines=show(p.output, p.top_m, progress=true)
    catch
        printstyled(p.output, "MODEL DISPLAY UNAVAILABLE\n", bold=true, color=:blue); println()
        return 3
    end
end

function model_obs_display(p)
    try
        println("Current MAP model")
        return lines=show(p.output, p.top_m, p.e, progress=true)
    catch
        printstyled(p.output, "MODEL DISPLAY UNAVAILABLE\n", bold=true, color=:blue); println()
        return 3
    end
end


text:
pref = Set()

f = open("hightemp.txt", "r")
for line in readlines(f)
  col1 = split(line, "\t")[1]
  push!(pref,col1)
end


text:
# m = (0,1,1,2)
# shuffle_24_bytes_to_32 = Val( ntuple(i->m[mod1(i,4)] +  3( (i-1)÷4 ),32) )
const shuffle_24_bytes_to_32 =  (
     0,  1,  1,  2,
     3,  4,  4,  5,
     6,  7,  7,  8,
     9, 10, 10, 11,
    12, 13, 13, 14,
    15, 16, 16, 17,
    18, 19, 19, 20,
    21, 22, 22, 23
)

# AVX 256-bit (32 bytes)
const shuffle_24_bytes_to_32_val = Val{shuffle_24_bytes_to_32}()
const shuffle_24_bytes_to_32_m4_val = Val{shuffle_24_bytes_to_32 .+ 4}()
const shuffle_24_bytes_to_32_m8_val = Val{shuffle_24_bytes_to_32 .+ 8}()
const choose_even_odd_ints_16 = Val( ntuple(i->isodd(i-1)*16 + i -1, 16) )

"""
    unpack_uint12_to_uint16(A::Vector{UInt8}, out::Vector{UInt16}, [ i ])

    Unpack 12-bit integers into 16-bit integers. Two 12-bit integers are packed consecutively into three bytes.
"""
@inline function unpack_uint12_to_uint16(
    A::SIMD.FastContiguousArray{UInt8,1},
    out::SIMD.FastContiguousArray{UInt16,1},
    i,
    offset = 0,
    shuffle_24_bytes_to_32::Val{T} = shuffle_24_bytes_to_32_val
    ) where T

    # Load 32 bytes (we only use the first 24)
    a = @inbounds vload(Vec{32,UInt8}, A, i-offset)

    a16 = reinterpret( Vec{16,UInt16}, shufflevector(a, shuffle_24_bytes_to_32) )

    a16 = shufflevector(a16 & 0xfff, a16 >> 4, choose_even_odd_ints_16)

    @inbounds vstore(a16, out, 1 + 2( ( i-1 ) ÷ 3) )
end
@inline unpack_uint12_to_uint16_offset_m4(A, out, i) = unpack_uint12_to_uint16(A, out, i, 4, shuffle_24_bytes_to_32_m4_val)
@inline unpack_uint12_to_uint16_offset_m8(A, out, i) = unpack_uint12_to_uint16(A, out, i, 8, shuffle_24_bytes_to_32_m8_val)

"""
    unpack_uint12_to_uint16(A::SIMD.FastContiguousArray{UInt8,1}, out::Vector{UInt16})

    Unpack entire array with a preallocated out buffer for compressed data
"""
function unpack_uint12_to_uint16(A::SIMD.FastContiguousArray{UInt8,1}, out::SIMD.FastContiguousArray{UInt16,1})
    unpack_uint12_to_uint16(A, out, 1)
    idx = 25:24:length(A)
    @inbounds for i = idx[1:end-1]
        unpack_uint12_to_uint16_offset_m4(A, out, i)
    end
    extra_bytes = rem(length(A), 3)
    unpack_uint12_to_uint16_offset_m8(A, out, length(A) - 23 - extra_bytes)
    if extra_bytes == 1
        out[length(A) * 2 ÷ 3] = A[end]
    elseif extra_bytes == 2
        out[length(A) * 2 ÷ 3] = reinterpret(UInt16, A[end-1:end])[1] & 0xfff
    end
    out
end

"""
    unpack_uint12_to_uint16(A::SIMD.FastContiguousArray{UInt8,1})

    Unpack entire array, allocates and returns output buffer
"""
function unpack_uint12_to_uint16(A::SIMD.FastContiguousArray{UInt8,1})
    out = Vector{UInt16}(undef, length(A) * 2 ÷ 3)
    unpack_uint12_to_uint16(A, out)
end

"""
    unpack_uint12_to_uint16(filepath::AbstractString)

    Load file at filepath, memory map it and unpack
"""
function unpack_uint12_to_uint16(filepath::AbstractString)
    ios = open(filepath, "r")
    A = Mmap.mmap(ios)
    unpack_uint12_to_uint16(A);
    close(ios)
end

function unpack_uint12_to_uint16_partitioned(inn, out = Array{UInt16,1}(undef, div(length(inn), 3, RoundUp)*2))
    partition_length = div( length(inn), Threads.nthreads(), RoundUp)
    # Ensure that partition length is a multiple of three
    partition_length += 3 - mod1(partition_length, 3)
    @debug "Partition Length: " partition_length partition_length ÷ 3 * 2
    inp = collect(Iterators.partition(inn, partition_length))
    outp = collect(Iterators.partition(out, partition_length ÷ 3 * 2))
    @inbounds Threads.@threads for i in eachindex(inp)
        unpack_uint12_to_uint16(inp[i], outp[i])
    end
    out
end

#=

=#
const shiftvec = reinterpret(Vec{16,UInt16},Vec{8,UInt32}(0x00040000))
const andvec = reinterpret(Vec{16,UInt16}, Vec{8,UInt32}(0xffff0fff))
@inline function unpack_uint12_to_uint16_one_shuffle(
    A::SIMD.FastContiguousArray{UInt8,1},
    out::SIMD.FastContiguousArray{UInt16,1},
    i,
    offset = 0,
    shuffle_24_bytes_to_32::Val{T} = shuffle_24_bytes_to_32_val,
    ) where T
    # Load 32 bytes (we only use the first 24)
    a = @inbounds vload(Vec{32,UInt8}, A, i-offset)

    a16 = reinterpret( Vec{16,UInt16}, shufflevector(a, shuffle_24_bytes_to_32) )

    a16 = a16 >> shiftvec & andvec

    @inbounds vstore(a16, out, 1 + 2( ( i-1 ) ÷ 3) )
end


text:
module LyceumMuJoCoViz

using Base: RefValue, @lock, @lock_nofail

# Stdlib
using Printf: @printf

# 3rd party
using GLFW: GLFW, Window, Key, Action, MouseButton, GetKey, RELEASE, PRESS, REPEAT
using PrettyTables: pretty_table
using BangBang: @set!!
using StaticArrays: SVector, MVector
using DocStringExtensions
using Observables: AbstractObservable, Observable, on, off
using FFMPEG

# Lyceum
using MuJoCo, MuJoCo.MJCore
using LyceumMuJoCo
import LyceumMuJoCo: reset!
using LyceumBase: LyceumBase, Maybe, AbsVec, AbsMat

export visualize

const FONTSCALE = MJCore.FONTSCALE_150 # can be 100, 150, 200
const MAXGEOM = 10000 # preallocated geom array in mjvScene
const MIN_REFRESHRATE = 30 # minimum rate when sim cannot run at the native refresh rate
const SIMGAMMA = 0.99
const RNDGAMMA = 0.9
const VIDFPS = 40

const RES_HD = (1280, 720)
const RES_FHD = (1920, 1080)
const RES_XGA = (1024, 768)
const RES_SXGA = (1280, 1024)

include("util.jl")
include("glfw.jl")
include("ratetimer.jl")
include("types.jl")
include("functions.jl")
include("modes.jl")
include("defaulthandlers.jl")

function __init__()
    if Threads.nthreads() == 1
        @warn "LyceumMuJoCoViz is designed to run multi-threaded, but the current Julia session was started with only one thread. Degraded performance will occur. To enable multi-threading, set JULIA_NUM_THREADS to a value greater than 1 before starting Julia."
    end
    return
end

"""
    $(TYPEDSIGNATURES)

Starts an interactive visualization of `model`, which can be a valid subtype of
`AbstractMuJoCoEnvironment`, a `MJSim` simulation, or the path to a valid MuJoCo model file.
The visualizer has several "modes" that allow you to visualize passive dynamics,
play back recorded trajectories, and run a controller interactively. The passive dynamics
mode depends only on `model` and is always available, while the other modes are specified
by the keyword arguments below.

# Keywords

- `trajectories`: a single trajectory or vector of trajectories, where each
    trajectory is an AbstractMatrix of states with size `(length(statespace(model)), T)` and
    `T` is the length of the trajectory. Note that each trajectory can have different length.
- `controller`: a callback function with the signature `controller(model)`, called at each
    timestep, that applies a control input to the system.
- `windowsize`: resolution of visualizer window. Defaults to an internal heuristic.

# Examples
```julia
using LyceumMuJoCo, LyceumMuJoCoViz
env = LyceumMuJoCo.HopperV2()
T = 100
states = Array(undef, statespace(env), T)
for t = 1:T
    step!(env)
    states[:, t] .= getstate(env)
end
visualize(
    env,
    trajectories = states,
    controller = env -> setaction!(env, rand(actionspace(env)))
)
```
"""
function visualize(
    model::Union{AbstractString,MJSim,AbstractMuJoCoEnvironment};
    trajectories::Union{Nothing,AbsMat,AbsVec{<:AbsMat}} = nothing,
    controller = nothing,
    windowsize::NTuple{2,Integer} = default_windowsize()
)
    model isa AbstractString && (model = MJSim(model))
    reset!(model)

    modes = EngineMode[PassiveDynamics()]
    !isnothing(trajectories) && push!(modes, Trajectory(trajectories))
    !isnothing(controller) && push!(modes, Controller(controller))

    run(Engine(windowsize, model, Tuple(modes)))
    return
end

function run(e::Engine)
    if e.phys.model isa AbstractMuJoCoEnvironment
        e.ui.reward = getreward(e.phys.model)
    end

    # render first frame before opening window
    prepare!(e)
    render(e)
    e.ui.refreshrate = GetRefreshRate()
    e.ui.lastrender = time()
    GLFW.ShowWindow(e.mngr.state.window)

    # run the simulation/mode in second thread
    modetask = Threads.@spawn runphysics(e)

    println(ASCII)
    println("Press \"F1\" to show the help message.")

    runui(e)
    wait(modetask)
    return
end

function runui(e::Engine)
    shouldexit = false
    trecord = 0.0
    try
        while !shouldexit
            @lock e.phys.lock begin
                GLFW.PollEvents()
                prepare!(e)
            end

            render(e)
            trender = time()

            rt = 1 / (trender - e.ui.lastrender)
            @lock e.ui.lock begin
                e.ui.refreshrate = RNDGAMMA * e.ui.refreshrate + (1 - RNDGAMMA) * rt
                e.ui.lastrender = trender
                shouldexit = e.ui.shouldexit |= GLFW.WindowShouldClose(e.mngr.state.window)
            end

            tnow = time()
            if e.ffmpeghandle !== nothing && tnow - trecord > 1 / e.min_refreshrate
                trecord = tnow
                recordframe(e)
            end

            yield()
        end
    finally
        @lock e.ui.lock begin
            e.ui.shouldexit = true
        end
        GLFW.DestroyWindow(e.mngr.state.window)
    end
    return
end

function prepare!(e::Engine)
    ui, p = e.ui, e.phys
    sim = getsim(p.model)
    p.model isa AbstractMuJoCoEnvironment && (ui.reward = getreward(p.model))
    mjv_updateScene(
        sim.m,
        sim.d,
        ui.vopt,
        p.pert,
        ui.cam,
        MJCore.mjCAT_ALL,
        ui.scn,
    )
    prepare!(ui, p, mode(e))
    return e
end

function render(e::Engine)
    w, h = GLFW.GetFramebufferSize(e.mngr.state.window)
    rect = mjrRect(Cint(0), Cint(0), Cint(w), Cint(h))
    mjr_render(rect, e.ui.scn, e.ui.con)
    e.ui.showinfo && overlay_info(rect, e)
    GLFW.SwapBuffers(e.mngr.state.window)
    return
end

function runphysics(e::Engine)
    p = e.phys
    ui = e.ui
    resettime!(p) # reset sim and world clocks to 0

    try
        while true
            shouldexit, lastrender, reversed, paused, refrate, = @lock_nofail ui.lock begin
                ui.shouldexit, ui.lastrender, ui.reversed, ui.paused, ui.refreshrate
            end

            if shouldexit
                break
            elseif (time() - lastrender) > 1 / e.min_refreshrate
                yield()
                continue
            else
                @lock p.lock begin
                    elapsedworld = time(p.timer)

                    # advance sim
                    if ui.paused
                        pausestep!(p, mode(e))
                    elseif ui.reversed && p.elapsedsim > elapsedworld
                        reversestep!(p, mode(e))
                        p.elapsedsim -= timestep(p.model)
                    elseif !ui.reversed && p.elapsedsim < elapsedworld
                        forwardstep!(p, mode(e))
                        p.elapsedsim += timestep(p.model)
                    end
                end
            end
        end
    finally
        @lock ui.lock begin
            ui.shouldexit = true
        end
    end

    return
end

end # module


text:
using Hecke

Base.:*(x::Hecke.NfRelNSElem{Nemo.nf_elem}) = x

function _get_simple_extension_and_maximal_order(K)
  @assert degree(base_ring(K)) == 1
  pol = K.pol
  k = length(pol)
  gensK = gens(K)
  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  basesofmaximalorders = Vector{elem_type(K)}[]
  discs = fmpz[]
  for i in 1:k
    f = pol[i]
    g = zero(Qx)
    for j in 1:f.length
      g += coeff(f.coeffs[j], 1)* x^Int(f.exps[i, j])
    end
    L, _ = number_field(g)
    OL = maximal_order(L)
    push!(discs, discriminant(OL))
    BOL = map(z -> z.elem_in_nf, basis(OL))
    push!(basesofmaximalorders, [ sum(coeff(b, j) * gensK[i]^j for j in 0:degree(L)) for b in BOL ])
  end
  bbb = vec([ prod(c) for c in Iterators.product(basesofmaximalorders...) ])
  Ksimpleabs, g = Hecke.absolute_simple_field(Ksimple)
  prime_divisors = Set{fmpz}()
  for i in 1:length(discs)
    for j in 1:(i - 1)
      ff = factor(gcd(discs[i], discs[j]))
      prime_divisors = union(prime_divisors, Set{fmpz}([ p for (p, e) in ff ]))
    end
  end
  OO = Order(Ksimpleabs, [ g\b for b in bbb ])
  for p in prime_divisors
    OO = pmaximal_overorder(OO, p)
  end
  OO.is_maximal = 1
  set_attribute!(Ksimpleabs, :maximal_order => OO)
  return Ksimpleabs
end

bounddisc = ARGS[1]
real = ARGS[2]
gtype = ARGS[3]
startcond = ARGS[4]
endcond = ARGS[5]

bounddisc = fmpz(eval(parse(bounddisc)))
gtype = convert(Vector{Int}, eval(parse(gtype)))
startcond = parse(Int, startcond)
endcond = parse(Int, endcond)
real = parse(Bool, real)

sprint_formatted(fmt, args...) = @eval @sprintf($fmt, $(args...))

Qx, x = PolynomialRing(QQ, "x")
K, a = NumberField(x - 1, "a")
O = maximal_order(K)

n=prod(gtype)
expo=lcm(gtype)

if length(ARGS) == 6
   total_cond = countlines(ARGS[6])
else
   conductors = sort!(Hecke.conductorsQQ(O, gtype, bounddisc))
   total_cond = length(conductors)
end

upper_index = min(total_cond, endcond)

l_conductors = Vector{Int}(upper_index - startcond + 1)

if length(ARGS) == 6
   print("Loading conductors ... ")
   open(ARGS[6]) do f
     k = 1
     for (i, lline) in enumerate(eachline(f))
       if i < startcond
         continue
       end
       if i > upper_index
         break
       end
       l_conductors[k] = parse(Int, lline)
       k = k + 1
     end
   end
   #conductors = readdlm(ARGS[6], Int)
   #conductors = conductors[:, 1]
   println("DONE")
end

width = length(string(total_cond))

#@show l_conductors

fields=Tuple{AnticNumberField, fmpz}[]
#  autos=Vector{NfRelNSToNfRelNSMor}[]

#Now, the big loop
for (i, k) in enumerate(l_conductors)
  print("Left: $(length(l_conductors) - i)\n")
  r,mr=Hecke.ray_class_groupQQ(O,k,!real,expo)
  if !Hecke._are_there_subs(r,gtype)
    continue
  end
  ls=subgroups(r,quotype=gtype, fun= (x, y) -> quo(x, y, false)[2])
  for s in ls
    C=Hecke.ray_class_field(mr, s)
    if Hecke._is_conductor_minQQ(C,n) && Hecke.discriminant_conductorQQ(O,C,k,bounddisc,n)
      L=number_field(C)
      LL = _get_simple_extension_and_maximal_order(L)
      LLdisc = discriminant(maximal_order(LL))

      if LLdisc != prod( p^e for (p, e) in C.absolute_discriminant)
        println("Ups")
        @show LLdisc
        @show prod( p^e for (p, e) in C.absolute_discriminant)
        println("==========================")
      end

      push!(fields, (Hecke.simplify(LL)[1], LLdisc))
    end
  end
end

if length(ARGS) == 6
file = ARGS[6] * "_" * sprint_formatted("%0$(width)d", startcond) * "-" * sprint_formatted("%0$(width)d", upper_index)
else
file = sprint_formatted("%0$(width)d", startcond) * "-" * sprint_formatted("%0$(width)d", upper_index)
end

@show length(fields)

Hecke._write_fields(fields, file)


text:
using AlgebraicPetri
using OrdinaryDiffEq
using Plots
using Catlab.Meta
using JSON

import OrdinaryDiffEq: ODEProblem
ODEProblem(p::LabelledReactionNet, t) = ODEProblem(vectorfield(p), concentrations(p), t, rates(p))

# help capture JSON of defined functions
macro capture(funcname, exname, ex)
    quote
        $(esc(exname)) = $(repr(strip_lines(ex, recurse=true)))
        $(esc(funcname)) = $ex
    end
end

@capture γ γ_text 1/14
@capture β β_text t->begin
    policy_days = [20,60,120] .+ 17
    contact_rate = 0.05
    pol = findfirst(x->t<=x, policy_days) # array of days when policy changes
    growth_rate = pol == 1 ? 0.0 : (2^(1/((pol-1)*5)) - 1) # growth rate depending on policy
    return (growth_rate + γ) / 990 * (1-contact_rate) # calculate rate of infection
end

sir_cset= LabelledReactionNet{Function, Float64}((:S=>990, :I=>10, :R=>0), (:inf, β)=>((:S, :I)=>(:I,:I)), (:rec, t->γ)=>(:I=>:R))

Graph(sir_cset)

prob = ODEProblem(sir_cset, (17.0, 120.0))
sol = OrdinaryDiffEq.solve(prob,Tsit5())
plot(sol)

## Getting Sharable JSON
sir_cset_string = LabelledReactionNet{String, Int}((:S=>990, :I=>10, :R=>0), (:inf, β_text)=>((:S, :I)=>(:I,:I)), (:rec, γ_text)=>(:I=>:R))
JSON.print(sir_cset_string.tables, 2)


text:
abstract type AbstractDataset end

Base.length(s::AbstractDataset) = size(s.values)[1]
get_data_resolution(s::AbstractDataset)::Dates.Millisecond = s.resolution
get_last_recorded_row(s::AbstractDataset) = s.last_recorded_row

"""
Return the timestamp from the data used in the last update
"""
get_update_timestamp(s::AbstractDataset) = s.update_timestamp

function set_last_recorded_row!(s::AbstractDataset, val::Int)
    s.last_recorded_row = val
    return
end

function set_update_timestamp!(s::AbstractDataset, val::Dates.DateTime)
    s.update_timestamp = val
    return
end

# Values field is accessed with dot syntax to avoid type instability

mutable struct DataFrameDataset <: AbstractDataset
    values::DataFrames.DataFrame
    # We use Array here to allow for overwrites when updating the state
    timestamps::Vector{Dates.DateTime}
    # Resolution is needed because AbstractDataset might have just one row
    resolution::Dates.Millisecond
    end_of_step_index::Int
    last_recorded_row::Int
    update_timestamp::Dates.DateTime
end

function DataFrameDataset(
    values::DataFrames.DataFrame,
    timestamps::Vector{Dates.DateTime},
    resolution::Dates.Millisecond,
    end_of_step_index::Int,
)
    return DataFrameDataset(
        values,
        timestamps,
        resolution,
        end_of_step_index,
        0,
        UNSET_INI_TIME,
    )
end

function DataFrameDataset(values::DataFrames.DataFrame)
    return DataFrameDataset(
        values,
        Vector{Dates.DateTime}(),
        Dates.Second(0.0),
        1,
        0,
        UNSET_INI_TIME,
    )
end

function make_system_state(
    values::DataFrames.DataFrame,
    timestamp::Dates.DateTime,
    resolution::Dates.Millisecond,
)
    return DataFrameDataset(values, [timestamp], resolution, 0, 1, UNSET_INI_TIME)
end

function get_dataset_value(s::DataFrameDataset, date::Dates.DateTime)
    s_index = find_timestamp_index(s.timestamps, date)
    if isnothing(s_index)
        error("Request time stamp $date not in the state")
    end
    return s.values[s_index, :]
end

function get_column_names(::OptimizationContainerKey, s::DataFrameDataset)
    return DataFrames.names(s.values)
end

function get_last_recorded_value(s::DataFrameDataset)
    if get_last_recorded_row(s) == 0
        error("The Dataset hasn't been written yet")
    end
    return s.values[get_last_recorded_row(s), :]
end

function get_end_of_step_timestamp(s::DataFrameDataset)
    return s.timestamps[s.end_of_step_index]
end

"""
Return the timestamp from most recent data row updated in the dataset. This value may not be the same as the result from `get_update_timestamp`
"""
function get_last_updated_timestamp(s::DataFrameDataset)
    last_recorded_row = get_last_recorded_row(s)
    if last_recorded_row == 0
        return UNSET_INI_TIME
    end
    return s.timestamps[last_recorded_row]
end

function get_value_timestamp(s::DataFrameDataset, date::Dates.DateTime)
    s_index = find_timestamp_index(s.timestamps, date)
    if isnothing(s_index)
        error("Request time stamp $date not in the state")
    end
    return s.timestamps[s_index]
end

function set_value!(s::DataFrameDataset, vals, index::Int)
    setindex!(s.values, vals, index, :)
    return
end

function set_value!(s::DataFrameDataset, vals::DataFrames.DataFrame, index::Int)
    @assert_op size(vals)[1] == 1
    set_value!(s, vals[1, :], index)
    return
end

# HDF5Dataset does not account of overwrites in the data. Values are written sequentially.
mutable struct HDF5Dataset <: AbstractDataset
    values::HDF5.Dataset
    column_dataset::HDF5.Dataset
    write_index::Int
    last_recorded_row::Int
    resolution::Dates.Millisecond
    initial_timestamp::Dates.DateTime
    update_timestamp::Dates.DateTime
end

HDF5Dataset(values, column_dataset, resolution, initial_time) =
    HDF5Dataset(values, column_dataset, 1, 0, resolution, initial_time, UNSET_INI_TIME)

function get_column_names(::OptimizationContainerKey, s::HDF5Dataset)
    return s.column_dataset[:]
end

"""
Return the timestamp from most recent data row updated in the dataset. This value may not be the same as the result from `get_update_timestamp`
"""
function get_last_updated_timestamp(s::HDF5Dataset)
    last_recorded_row = get_last_recorded_row(s)
    if last_recorded_row == 0
        return UNSET_INI_TIME
    end
    return s.initial_timestamp + s.resolution * (last_recorded_row - 1)
end

function get_value_timestamp(s::HDF5Dataset, date::Dates.DateTime)
    s_index = find_timestamp_index(s.timestamps, date)
    if isnothing(s_index)
        error("Request time stamp $date not in the state")
    end
    return s.initial_timestamp + s.resolution * (s_index - 1)
end

function set_value!(s::HDF5Dataset, vals, index::Int)
    # Temporary while there is no implementation of caching of em_data
    _write_dataset!(s.values, vals, index:index)
    return
end


text:
using HuggingFaceApi
using Documenter

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

makedocs(;
    modules=[HuggingFaceApi],
    authors="chengchingwen <adgjl5645@hotmail.com> and contributors",
    repo="https://github.com/chengchingwen/HuggingFaceApi.jl/blob/{commit}{path}#{line}",
    sitename="HuggingFaceApi.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://chengchingwen.github.io/HuggingFaceApi.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/chengchingwen/HuggingFaceApi.jl",
)


text:
using IntervalArithmetic
using DifferentialEquations
using LinearAlgebra
using TaylorIntegration
using DiffEqBase
# setprecision(100)
function Flo(du,u,p,t)
    h=@interval(.125)
    A11=(-1) *(1+4. *h^2+4. *h*cos(2. *t)).^(-1/2).*sin(2. *t)
    A12=(1+4. *h^2+(1+4. *h^2+4. *h*cos(2. *t)).^(1/2)+2. *h*cos(2. *t).*(2+(1+4. *h^2+ 4. *h*cos(2. *t)).^(1/2))).^(-1) *((-1) *cos(2. *t).*(1+(-8).*h^2+(1+4. *h^2+4. *h*cos(2. *t)).^(1/2))+h*(1+(-4).*h+16. *h^2+2. *(1+ 4. *h^2+4. *h*cos(2. *t)).^(1/2)+(-8).*h*(1+4. *h^2+4. *h*cos(2. * t)).^(1/2)+(-1) *cos(4. *t)))
    A21=(1+4. *h^2+(1+4. *h^2+4. *h*cos(2. *t) ).^(1/2)+2. *h*cos(2. *t).*(2+(1+4. *h^2+4. *h*cos(2. *t)).^(1/2))).^(-1) *(h+4. *h^2. *(1+4. *h)+2. *h*(1+4. *h).*(1+4. *h^2+4. *h*cos(2. *t)).^(1/2)+(-1) *cos(2. *t).*(1+(-8).*h^2+(1+4. *h^2+4. *h*cos( 2. *t)).^(1/2))+(-1) *h*cos(4. *t))
    A22=(1+4. *h^2+4. *h*cos(2. *t)).^( .-1/2).*sin(2. *t)

    du[1]=A11*u[1]+A12*u[2]
    du[2]=A21*u[1]+A22*u[2]
end
#

T_end=pi
tspan=(0.,T_end)
tol=1e-16

u0=[1.,0.]
prob = ODEProblem(Flo,u0,tspan)

# t,u=taylorinteg(Flo,u0, tspan[1],tspan[2],25,1e-50, maxsteps=1_000_000_000)
t,u=solve(prob,TaylorMethod(10),abstol=tol,maxiters=1_000_000_000)
diff=sqrt((u[1,end]-u0[1])^2+(u[2,end]-u0[2])^2)
# diff=sqrt((u[1,end]-u1[1,1])^2+(u1[2,end]-u1[2,1])^2)
# println(diff)
# using IntervalArithmetic
# using DifferentialEquations
# using LinearAlgebra
# using TaylorIntegration
# using DiffEqBase
# setprecision(100)
# @taylorize function Flo(du,u,p,t)
#     h=@interval(.125)
#     A11=(-1) *(1+4. *h^2+4. *h*cos(2. *t)).^(-1/2).*sin(2. *t)
#     A12=(1+4. *h^2+(1+4. *h^2+4. *h*cos(2. *t)).^(1/2)+2. *h*cos(2. *t).*(2+(1+4. *h^2+ 4. *h*cos(2. *t)).^(1/2))).^(-1) *((-1) *cos(2. *t).*(1+(-8).*h^2+(1+4. *h^2+4. *h*cos(2. *t)).^(1/2))+h*(1+(-4).*h+16. *h^2+2. *(1+ 4. *h^2+4. *h*cos(2. *t)).^(1/2)+(-8).*h*(1+4. *h^2+4. *h*cos(2. * t)).^(1/2)+(-1) *cos(4. *t)))
#     A21=(1+4. *h^2+(1+4. *h^2+4. *h*cos(2. *t) ).^(1/2)+2. *h*cos(2. *t).*(2+(1+4. *h^2+4. *h*cos(2. *t)).^(1/2))).^(-1) *(h+4. *h^2. *(1+4. *h)+2. *h*(1+4. *h).*(1+4. *h^2+4. *h*cos(2. *t)).^(1/2)+(-1) *cos(2. *t).*(1+(-8).*h^2+(1+4. *h^2+4. *h*cos( 2. *t)).^(1/2))+(-1) *h*cos(4. *t))
#     A22=(1+4. *h^2+4. *h*cos(2. *t)).^( .-1/2).*sin(2. *t)
#
#     du[1]=A11*u[1]+A12*u[2]
#     du[2]=A21*u[1]+A22*u[2]
# end
# #
#
# T_end=BigFloat(pi)
# tspan=(BigFloat("0"),T_end)
# tol=BigFloat("1e-25")
#
# u0=[BigFloat("1"),BigFloat("0")]
# prob = ODEProblem(Flo,u0,tspan)
#
# # t,u=taylorinteg(Flo,u0, tspan[1],tspan[2],25,1e-50, maxsteps=1_000_000_000)
# t,u=solve(prob,TaylorMethod(25),abstol=tol,maxiters=1_000_000_000_000)
# diff=sqrt((u[1,end]-u0[1])^2+(u[2,end]-u0[2])^2)
# # diff=sqrt((u[1,end]-u1[1,1])^2+(u1[2,end]-u1[2,1])^2)
# println(diff)


text:
include("LUdecomp.jl")
include("Lsolve.jl")
include("Usolve.jl")

#we write a function that solves Ax = b for x by LU decomposition
function Asolve(A, b)
	L, U = LUdecomp(A)
	y = Lsolve(L, b)
	x = Usolve(U, y)
	return x
end

max_stars_repo_path: 1a/LU Decomposition/Asolve.jl

max_issues_repo_path: 1a/LU Decomposition/Asolve.jl

max_forks_repo_path: 1a/LU Decomposition/Asolve.jl


text:
using FancyPictures
using Test

@testset "FancyPictures.jl" begin
    # Write your own tests here.
end


text:
# src/qccdevctrlutils.jl
# Created by Alejandro Villoria and Anabel Ovide, 21 May, 2021
# MIT license
# Sub-module QCCDevCtrl

module QCCDev_Utils

export giveZone

using ..QCCDevControl_Types

"""
Given a zone ID, return the zone from the QCCDevControl object.
Returns `nothing` if no zone with that ID is found.
"""
function giveZone(qdc ::QCCDevControl,
    id ::Symbol)::Union{GateZone, AuxZone, LoadingZone, Junction, Nothing}

    zone = get(qdc.gateZones, id, nothing)
    if isnothing(zone)
        zone = get(qdc.junctions, id, nothing)
    end
    if isnothing(zone)
        zone = get(qdc.auxZones, id, nothing)
    end
    if isnothing(zone)
        zone = get(qdc.loadingZones, id, nothing)
    end
    return zone
end

end # Module

max_stars_repo_name: TartuQC/trapped-ion-qccd-simulator

max_issues_repo_name: TartuQC/trapped-ion-qccd-simulator

max_forks_repo_name: TartuQC/trapped-ion-qccd-simulator


text:
for (elty, relty, ext) in ((:ElInt, :ElInt, :i),
                           (:Float32, :Float32, :s),
                           (:Float64, :Float64, :d),
                           (:Complex64, :Float32, :c),
                           (:Complex128, :Float64, :z))

    for (mat, sym) in ((:Matrix, "_"),
                       (:DistMatrix, "Dist_"),
                       (:DistMultiVec, "DistMultiVec_"))
        @eval begin
            # Bernoulli
            function bernoulli!(A::$mat{$elty}, m::Integer = size(A, 1), n::Integer = 1, p::Real = 0.5)
                err = ccall(($(string("ElBernoulli", sym, ext)), libEl), Cuint,
                    (Ptr{Void}, ElInt, ElInt, Float64),
                    A.obj, m, n, Float64(p))
                err == 0 || throw(ElError(err))
                return A
            end

            # Gaussian
            function gaussian!(A::$mat{$elty}, m::Integer = size(A, 1), n::Integer = 1,
                               mean::Number = 0, stddev::Number = 1)
                err = ccall(($(string("ElGaussian", sym, ext)), libEl), Cuint,
                    (Ptr{Void}, ElInt, ElInt, $elty, $relty),
                    A.obj, m, n, mean, stddev)
                err == 0 || throw(ElError(err))
                return A
            end

            # Ones
            function ones!(A::$mat{$elty}, m::Integer = size(A, 1), n::Integer = 1)
                err = ccall(($(string("ElOnes", sym, ext)), libEl), Cuint,
                    (Ptr{Void}, ElInt, ElInt),
                    A.obj, m, n)
                err == 0 || throw(ElError(err))
                return A
            end

            # Uniform
            function uniform!(A::$mat{$elty}, m::Integer = size(A, 1), n::Integer = 1,
                              center::Number = 0, radius::Number = 1)
                err = ccall(($(string("ElUniform", sym, ext)), libEl), Cuint,
                    (Ptr{Void}, ElInt, ElInt, $elty, $relty),
                    A.obj, m, n, center, radius)
                err == 0 || throw(ElError(err))
                return A
            end

            # Zeros
            function zeros!(A::$mat{$elty}, m::Integer = size(A, 1), n::Integer = 1)
                err = ccall(($(string("ElZeros", sym, ext)), libEl), Cuint,
                    (Ptr{Void}, ElInt, ElInt),
                    A.obj, m, n)
                err == 0 || throw(ElError(err))
                return A
            end
        end

        if elty == :Complex64 || elty == :Complex128
            # Uniform
            @eval begin
                function foxLi!(A::$mat{$elty}, n::Integer = size(A, 1), ω::Real = 1.0)
                    err = ccall(($(string("ElFoxLi", sym, ext)), libEl), Cuint,
                        (Ptr{Void}, ElInt, $relty),
                        A.obj, n, ω)
                    err == 0 || throw(ElError(err))
                    return A
                end
            end
        end
    end
end

for (elty, relty, ext) in ((:Float32, :Float32, :s),
                           (:Float64, :Float64, :d),
                           (:Complex64, :Float32, :c),
                           (:Complex128, :Float64, :z))

    for (mat, sym) in ((:Matrix, "_"),
                       (:DistMatrix, "Dist_"),
                       (:SparseMatrix, "Sparse_"),
                       (:DistSparseMatrix, "DistSparse_"))
        @eval begin
            # Helmholtz
            function helmholtz!(A::$mat{$elty}, nx::Integer; shift::Number = 0)
                err = ccall(($(string("ElHelmholtz1D", sym, ext)), libEl), Cuint,
                    (Ptr{Void}, ElInt, $elty),
                    A.obj, nx, shift)
                err == 0 || throw(ElError(err))
                return A
            end

            function helmholtz!(A::$mat{$elty}, nx::Integer, ny::Integer; shift::Number = 0)
                err = ccall(($(string("ElHelmholtz2D", sym, ext)), libEl), Cuint,
                    (Ptr{Void}, ElInt, ElInt, $elty),
                    A.obj, nx, ny, shift)
                err == 0 || throw(ElError(err))
                return A
            end

            function helmholtz!(A::$mat{$elty}, nx::Integer, ny::Integer, nz::Integer; shift::Number = 0)
                err = ccall(($(string("ElHelmholtz3D", sym, ext)), libEl), Cuint,
                    (Ptr{Void}, ElInt, ElInt, ElInt, $elty),
                    A.obj, nx, ny, nz, shift)
                err == 0 || throw(ElError(err))
                return A
            end
        end
    end
end


text:
const two27 = Float64(1<<27)
const two12 = Float32(1<<12)
const two6  = Float16(1<<6)

"""
    extractscalar

Split a value relative to an integral power of 2

hi is the high order part of p, lo is the residual part
hi + lo === p
"""

for (T,P) in ((:Float64, 53), (:Float32, 24), (:Float16, 11))
  @eval begin
    @inline function extractscalar(value::$T)
        powoftwo =  ldexp(one($T), ($P- exponent(value)) >> 1) # P == precision(Float64)
        return extractscalar(powoftwo, value)
    end
  end
end

function extractscalar(powoftwo::T, value::T) where {T}
   hi = (powoftwo + value) - powoftwo
   lo = value - hi
   return hi, lo
end

"""
    two_cube(a)

Computes `s = fl(a*a*a)` and `e1 = err(a*a*a)`.
"""
function two_cube(value::T) where {T}
    hi, lo = extractscalar(value)
    hhh  = three_prod(hi, hi, hi)
    hhl  = three_prod(hi, hi, 3*lo)
    hll  = three_prod(hi, 3*lo, lo)
    lll  = three_prod(lo, lo, lo)

    himh = four_sum(hhh[1], hhl[1], hhh[2], hll[1])
    mllo = four_sum(hhl[2], hll[2], lll[1], lll[2])
    hilo = four_sum(himh[1], himh[2], mllo[1], mllo[2])

    return hilo[1], hilo[2]
end

"""
    three_cube(a)

Computes `s = fl(a*a*a)` and `e1 = err(a*a*a), e2 = err(e1)`.
"""
function three_cube(value::T) where {T}
    hi, lo = extractscalar(value)
    hhh  = three_prod(hi, hi, hi)
    hhl  = three_prod(hi, hi, 3*lo)
    hll  = three_prod(hi, 3*lo, lo)
    lll  = three_prod(lo, lo, lo)

    himh = four_sum(hhh[1], hhl[1], hhh[2], hll[1])
    mllo = four_sum(hhl[2], hll[2], lll[1], lll[2])
    hilo = four_sum(himh[1], himh[2], mllo[1], mllo[2])

    return hilo[1], hilo[2], hilo[3]
end

"""
    four_cube(a)

Computes `s = fl(a*a*a)` and `e1 = err(a*a*a), e2 = err(e1), e3 = err(e2)`.
"""
function three_cube(value::T) where {T}
    hi, lo = extractscalar(value)
    hhh  = three_prod(hi, hi, hi)
    hhl  = three_prod(hi, hi, 3*lo)
    hll  = three_prod(hi, 3*lo, lo)
    lll  = three_prod(lo, lo, lo)

    himh = four_sum(hhh[1], hhl[1], hhh[2], hll[1])
    mllo = four_sum(hhl[2], hll[2], lll[1], lll[2])
    hilo = four_sum(himh[1], himh[2], mllo[1], mllo[2])

    return hilo
end


text:
include("common.jl")
bcfl =  TwoAgents.BCFL21()

for i in workers()
    remotecall_fetch(i, include, "common.jl")
end

TwoAgents.brutal_solution(bcfl; maxiter=1000)

max_stars_repo_name: NYUEcon/GrowthModels

max_issues_repo_name: NYUEcon/GrowthModels

max_forks_repo_name: NYUEcon/GrowthModels


text:
using ViewerGL
GL = ViewerGL

points = rand(1000,3)
GL.VIEW([
      GL.GLFrame,
      GL.GLPoints(points)
]);

points = rand(1000,2)
GL.VIEW([
      GL.GLFrame2,
      GL.GLPoints(points)
]);


text:
	#
	#	Script for calculation of far-field pattern in the Dirichlet case (prolate)
	#

	# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	# User configurable parameters
	# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

	# Physical parameters
	k_0 = 4 ; # Wave number in media 0
	a = 1 ; # Major semiaxis spheroid
	b = 0.8 ; # Minor semiaxis spheroid
	theta_inc = pi/4 ; # Incidence angle in radians

	# Software parameters
	M = 8 ;
	method = 2 ;
	delta_eta = 0.005 ;

	# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	# Calculation
	# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	# Conversion theta to eta
	eta_inc = cos( theta_inc ) ;
	# Coefficient evaluation
	Amn = Coeff_ProDN( k_0, a, b, M, method, eta_inc, "D" ) ;

	# Far-Field pattern calculation
	Pattern = Pattern_LiquidPro( k_0, a, b, M, method, eta_inc, Amn, delta_eta ) ;

	# Saving to disk
	writedlm( "Out.Pattern.dat", Pattern , '\t' ) ;


text:
using Omega
using Flux

function testsum(;kwargs...)
  x = logistic(1.0, 1.0, (10,))
  y = sum(x)
  ΩT = Omega.SimpleΩ{Vector{Int}, Flux.TrackedArray{Float64, 1, Array{Float64,1}}}
  samples = rand(y, y ==ₛ -1.0; alg = HMCFAST, ΩT = ΩT, kwargs...)
end

testsum()


text:
#!/usr/bin/env julia

# Simple graphs
###############

run(pipeline(`dot -Tjson0 graphviz_digraph.dot`, stdout="graphviz_digraph.json"))
run(pipeline(`neato -Tjson0 graphviz_graph.dot`, stdout="graphviz_graph.json"))

# Wiring diagrams
#################

using Catlab.Graphics
using Catlab.Graphics.Graphviz: run_graphviz

diagram = include("graphviz_wiring_diagram.jl")
graph = to_graphviz(diagram)
open("graphviz_wiring_diagram.json", "w") do io
  run_graphviz(io, graph, format="json0")
end


text:
@testset "Trajectory" begin

    s = get_simulation("flamingo", "flat_2D_lc", "flat")
    ref_traj = deepcopy(ContactImplicitMPC.get_trajectory(s.model, s.env,
        joinpath(module_dir(), "src/dynamics/flamingo/gaits/gait_forward_36_4.jld2"),
        load_type = :split_traj_alt))

    N_repeat = 5
    N_sample = 1
    dupl_traj = repeat_ref_traj(ref_traj, N_repeat, idx_shift=[1])
    @test all(ContactImplicitMPC.tracking_error(ref_traj, dupl_traj, N_sample, idx_shift=[1]) .== 0.0)

end


text:
using Documenter
using Luna

modulesdir = joinpath(Luna.Utils.lunadir(), "docs", "src", "modules")

makedocs(
    sitename = "Luna Documentation",
    authors = "Christian Brahms and John C. Travers",
    pages = Any[
        "Home" => "index.md",
        "The numerical model" => [
            "General description" => "model/model.md",
            "Modal decompositions" => "model/modal_decompositions.md",
            "Nonlinear responses" => "model/nonlinear_responses.md"
        ],
        "The simple interface" => "interface.md",
        "Parameter scans" => "scans.md",
        "Modules" => [
            "$(split(fi, ".")[1]).jl" => "modules/$fi" for fi in readdir(modulesdir)
        ],
    ],
    format = Documenter.HTML(
        prettyurls = false
    )
)

deploydocs(
    repo = "github.com/LupoLab/Luna.jl.git",
)


text:
module TestLazyImports

using LazyImports
using Test

module TestShim
    using LazyImports
    @shim_import Test="8dfed614-e22c-5e08-85e1-65c5234f0b40"
end

module TestLazy
    using LazyImports
    @lazy_import Test="8dfed614-e22c-5e08-85e1-65c5234f0b40"
end

@testset "LazyImports.jl" begin
    @test TestShim.Test.detect_ambiguities === Test.detect_ambiguities
    @test TestLazy.Test.detect_ambiguities === Test.detect_ambiguities
end

end  # module


text:
function ConstantVelocity(; name, v = 1.0)
    @named oneport = OnePort()
    @unpack u = oneport
    ps = @parameters x=x # TODO: does this cause naming issue?
    eqs = [
        v ~ u
    ]
    extend(ODESystem(eqs, t, [], ps; name), oneport)
end

function Force(; name, F = 1.0)
    @named oneport = OnePort()
    @unpack u, f = oneport
    ps = @parameters F=F
    eqs = [
        F ~ f
    ]
    extend(ODESystem(eqs, t, [], ps; name), oneport)
end


text:
# include("AngularMomentum.jl")
# include("Pseudopotentials.jl")

using Combinatorics

abstract type HaldaneSphereStatesList end

struct HaldaneSphereStatesListLLL <: HaldaneSphereStatesList
    N::Int
    Nphi::Int
    Lz::AngMom

    possible_mz::Vector{Float64}

    orblist::Vector{Vector{Int}}
    nlist::Vector{Vector{Int}} # FIXME: use BitVector
    Ilist::Vector{Int}

    function HaldaneSphereStatesListLLL(N::Int, Nphi::Int, Lz::AngMom)
        Norb = Nphi + 1

        @assert N <= Norb
        @assert isangmom(Lz)

        Q = Nphi/2
        possible_mz = -Q:Q

        println("Building states list for N = $(N), Nphi = $(Nphi), Lz = $(Lz) ...")
        tic()

        orblist = Vector{Int}[] # occupied orbitals (ordered set containing N elements from 1:Norb)
        nlist = Vector{Int}[] # bit string of orbital occupations
        Ilist = Int[] # integer representative of state
        for orbs in combinations(1:Norb, N)
            if sum(possible_mz[orbs]) == Lz
                # orblist
                @assert issorted(orbs)
                push!(orblist, orbs)
                # nlist
                ns = nvec(orbs, Norb)
                push!(nlist, ns)
                # Ilist
                push!(Ilist, I(ns))
            end
        end

        sort_inds = sortperm(Ilist)
        orblist = orblist[sort_inds]
        nlist = nlist[sort_inds]
        Ilist = Ilist[sort_inds]

        println("Finished building states list of size $(length(Ilist))")
        toc()
        println()

        return new(N, Nphi, Lz, possible_mz, orblist, nlist, Ilist)
    end
end

get_N(states::HaldaneSphereStatesListLLL) = states.N
get_Nphi(states::HaldaneSphereStatesListLLL) = states.Nphi
get_Norb(states::HaldaneSphereStatesListLLL) = states.Nphi + 1
get_ell(states::HaldaneSphereStatesListLLL) = states.Nphi/2 # == angular momentum in the effective LLL problem == monopole strength Q in LLL

dim(states::HaldaneSphereStatesList) = length(states.Ilist)

mzvec(states::HaldaneSphereStatesListLLL, j::Int) = states.possible_mz[states.orblist[j]]

function nvec(orbs::Vector{Int}, Norb::Int)
    ns = zeros(Int, Norb)
    ns[orbs] = 1
    return ns
end

I(ns::Vector{Int}) = sum([ns[alpha]*2^(alpha-1) for alpha in 1:length(ns)])

ph_conjugate(ns::Vector{Int}) = 1 - ns

function ph_conjugate{T}(psi::Vector{T}, states::HaldaneSphereStatesListLLL)
    @assert length(psi) == dim(states)
    psi_ph = zeros(T, length(psi))
    for i in 1:length(psi)
        psi_ph[i] = psi[searchsorted(states.Ilist, I(ph_conjugate(states.nlist[i])))[1]]
    end
    return psi_ph
end

function build_Lsquared(states::HaldaneSphereStatesListLLL)
    println("Building L^2 matrix of size $(dim(states)) x $(dim(states)) ...")
    tic()

    M = get_Norb(states)
    Q = get_ell(states)

    rows = Int[]
    cols = Int[]
    LpLmvals = Float64[]
    LmLpvals = Float64[]
    @inbounds begin
    for j in 1:dim(states)
        if j == 1 || mod(j, ceil(0.05dim(states))) == 0
            println("Working on LpLm, column j = $(j): $(j)/$(dim(states)) = $(100j/dim(states))%")
        end
        ns = states.nlist[j]
        # diagonal terms of LpLm and LmLp
        push!(rows, j)
        push!(cols, j)
        push!(LpLmvals, sum([Aplus(Q, m(alpha, Q) - 1) * Aminus(Q, m(alpha, Q)) * ns[alpha] for alpha in 1:M])) # FIXME: don't use list comprehension
        push!(LmLpvals, sum([Aminus(Q, m(alpha, Q) + 1) * Aplus(Q, m(alpha, Q)) * ns[alpha] for alpha in 1:M])) # FIXME: don't use list comprehension
        # off-diagonal terms of LpLm and LmLp (they are identical)
        for alpha1 in 1:M-1
            for alpha2 in 2:M
                if (alpha1 != alpha2) && (alpha1 + 1 != alpha2 - 1)
                    if ns[alpha1] == 1 && ns[alpha2] == 1
                        sgn = alpha1 < alpha2 ? 1 : -1
                        ns_new = copy(ns)
                        ns_new[alpha1] = 0
                        ns_new[alpha2] = 0
                        if ns_new[alpha1 + 1] == 0 && ns_new[alpha2 - 1] == 0
                            sgn *= alpha1 + 1 < alpha2 - 1 ? 1 : -1
                            val = Aplus(Q, m(alpha1, Q)) * Aminus(Q, m(alpha2, Q))
                            if val != 0
                                p1 = sum(ns[1:alpha1-1])
                                p2 = sum(ns[1:alpha2-1])
                                p1p1 = sum(ns_new[1:(alpha1+1)-1])
                                p2m1 = sum(ns_new[1:(alpha2-1)-1])
                                sgn *= (-1)^(p1 + p2 + p1p1 + p2m1) # sign of matrix element from anticommutation relations
                                ns_new[alpha1 + 1] = 1
                                ns_new[alpha2 - 1] = 1
                                i = searchsorted(states.Ilist, I(ns_new))[1] # binary search to figure out which row to populate
                                push!(rows, i)
                                push!(cols, j)
                                push!(LpLmvals, -sgn * val)
                                push!(LmLpvals, -sgn * val)
                            end
                        end
                    end
                end
            end
        end
    end
    end
    LpLm = sparse(rows, cols, LpLmvals, dim(states), dim(states))
    LmLp = sparse(rows, cols, LmLpvals, dim(states), dim(states))

    cols = Int[]
    Lz2vals = Float64[]
    for j in 1:dim(states)
        if j == 1 || mod(j, ceil(0.10dim(states))) == 0
            println("Working on Lz^2, column j = $(j): $(j)/$(dim(states)) = $(100j/dim(states))%")
        end
        val = float(sum([m(alpha, Q) * states.nlist[j][alpha] for alpha in 1:M]))^2
        if val != 0
            push!(cols, j)
            push!(Lz2vals, val)
        end
    end
    Lz2 = sparse(cols, cols, Lz2vals, dim(states), dim(states))

    L2 = 0.5 * (LpLm + LmLp) + Lz2

    toc()
    println()

    return L2
end

abstract type HaldaneSphereSetup end

struct HaldaneSphereSetupLLL <: HaldaneSphereSetup
    states::HaldaneSphereStatesListLLL
    Lsquared::SparseMatrixCSC{Float64, Int64}

    function HaldaneSphereSetupLLL(N::Int, Nphi::Int, Lz::AngMom)
        # FIXME: add some asserts

        states = HaldaneSphereStatesListLLL(N, Nphi, Lz)
        Lsquared = build_Lsquared(states)

        return new(states, Lsquared)
    end
end

get_N(setup::HaldaneSphereSetupLLL) = get_N(setup.states)
get_Nphi(setup::HaldaneSphereSetupLLL) = get_Nphi(setup.states)
get_Norb(setup::HaldaneSphereSetupLLL) = get_Norb(setup.states)
get_ell(setup::HaldaneSphereSetupLLL) = get_ell(setup.states)

dim(setup::HaldaneSphereSetupLLL) = dim(setup.states)

function eig_Lsquared(Lsquared::SparseMatrixCSC{Float64, Int64})
    @assert ishermitian(Lsquared)
    println("Diagonalizing L^2 matrix of size $(size(Lsquared)[1]) x $(size(Lsquared)[2]) with eig ...")
    tic()
    evals, evecs = eig(Hermitian(full(Lsquared)))
    toc()
    println()
    return Lqn(evals), evecs
end

eig_Lsquared(setup::HaldaneSphereSetupLLL) = eig_Lsquared(setup.Lsquared)

# eigs on Lsquared is sketchy; use with caution..
function eigs_Lsquared(Lsquared::SparseMatrixCSC{Float64, Int64}, nev::Int)
    @assert ishermitian(Lsquared)
    println("Diagonalizing L^2 matrix of size $(size(Lsquared)[1]) x $(size(Lsquared)[2]) with eigs, nev = $(nev) ...")
    tic()
    evals, evecs = eigs(Lsquared, which=:SR, nev=nev)
    toc()
    println()
    return Lqn(evals), evecs
end

eigs_Lsquared(setup::HaldaneSphereSetupLLL, nev::Int) = eigs_Lsquared(setup.Lsquared, nev)

abstract type HaldaneSphereHami end

mutable struct HaldaneSphereHamiLLL <: HaldaneSphereHami
    setup::HaldaneSphereSetupLLL

    Vmat::Array{Float64, 4}
    H::SparseMatrixCSC{Float64, Int}

    energies::Vector{Float64}
    eigenstates::Matrix{Float64}

    spectrum::Matrix{Float64} # [L, E] pairs as rows in a 2-column matrix

    # use Vmat as input (most general case)
    function HaldaneSphereHamiLLL(setup::HaldaneSphereSetupLLL, Vmat::Array{Float64, 4})
        Norb = get_Norb(setup)
        @assert size(Vmat)[1] == size(Vmat)[2] == size(Vmat)[3] == size(Vmat)[4]
        @assert size(Vmat)[1] == Norb

        println("Building Hami matrix of size $(dim(setup)) x $(dim(setup)) ...")
        tic()
        rows = Int[]
        cols = Int[]
        Hvals = Float64[]
        @inbounds begin
        for j in 1:dim(setup)
            if j == 1 || mod(j, ceil(0.05dim(setup))) == 0
                println("Working on H, column j = $(j): $(j)/$(dim(setup)) = $(100j/dim(setup))%")
            end
            ns = setup.states.nlist[j]
            for alpha2 in 1:Norb
                for alpha1 in 1:alpha2-1
                    for alpha3 in 1:Norb
                        for alpha4 in 1:alpha3-1
                            if alpha1 + alpha2 == alpha3 + alpha4 # Lz conservation [mz = alpha - (Q + 1) = alpha + const.]
                                @assert alpha1 < alpha2 && alpha4 < alpha3
                                if ns[alpha3] == 1 && ns[alpha4] == 1
                                    ns_new = copy(ns)
                                    ns_new[alpha3] = 0
                                    ns_new[alpha4] = 0
                                    if ns_new[alpha1] == 0 && ns_new[alpha2] == 0
                                        V1234 = Vmat[alpha1, alpha2, alpha3, alpha4]
                                        V1243 = Vmat[alpha1, alpha2, alpha4, alpha3]
                                        V2134 = Vmat[alpha2, alpha1, alpha3, alpha4]
                                        V2143 = Vmat[alpha2, alpha1, alpha4, alpha3]
                                        val = V1234 - V1243 - V2134 + V2143
                                        if val != 0
                                            p1 = sum(ns_new[1:alpha1-1])
                                            p2 = sum(ns_new[1:alpha2-1])
                                            p3 = sum(ns[1:alpha3-1])
                                            p4 = sum(ns[1:alpha4-1])
                                            sgn = (-1)^(p1 + p2 + p3 + p4) # sign of matrix element from anticommutation relations
                                            ns_new[alpha1] = 1
                                            ns_new[alpha2] = 1
                                            i = searchsorted(setup.states.Ilist, I(ns_new))[1] # binary search to figure out which row to populate
                                            push!(rows, i)
                                            push!(cols, j)
                                            push!(Hvals, sgn * 0.5val)
                                        end
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
        end
        H = sparse(rows, cols, Hvals, dim(setup), dim(setup))
        toc()
        println()

        return new(setup, Vmat, H)
    end
end

# compute Clebsch–Gordan coefficients from scratch
function HaldaneSphereHamiLLL(setup::HaldaneSphereSetupLLL, pp::SphericalPseudopotentials)
    return HaldaneSphereHamiLLL(setup, build_Vmat(pp, build_CG_table(get_ell(setup))))
end

# precomputed Clebsch–Gordan coefficients as input
function HaldaneSphereHamiLLL(setup::HaldaneSphereSetupLLL, pp::SphericalPseudopotentials, CG_table::Array{Vector{Float64}, 3})
    return HaldaneSphereHamiLLL(setup, build_Vmat(pp, CG_table))
end

# compute Clebsch–Gordan coefficients from scratch
function HaldaneSphereHamiLLL(N::Int, Nphi::Int, Lz::AngMom, pp::SphericalPseudopotentials)
    return HaldaneSphereHamiLLL(HaldaneSphereSetupLLL(N, Nphi, Lz), build_Vmat(pp, build_CG_table(get_ell(setup))))
end

# precomputed Clebsch–Gordan coefficients as input
function HaldaneSphereHamiLLL(N::Int, Nphi::Int, Lz::AngMom, pp::SphericalPseudopotentials, CG_table::Array{Vector{Float64}, 3})
    return HaldaneSphereHamiLLL(HaldaneSphereSetupLLL(N, Nphi, Lz), build_Vmat(pp, CG_table))
end

get_N(Hami::HaldaneSphereHamiLLL) = get_N(Hami.setup)
get_Nphi(Hami::HaldaneSphereHamiLLL) = get_Nphi(Hami.setup)
get_Norb(Hami::HaldaneSphereHamiLLL) = get_Norb(Hami.setup)
get_ell(Hami::HaldaneSphereHamiLLL) = get_ell(Hami.setup)

dim(Hami::HaldaneSphereHamiLLL) = dim(Hami.setup)

function eig!(Hami::HaldaneSphereHamiLLL)
    @assert ishermitian(Hami.H)
    println("Diagonalizing Hami of size $(dim(Hami)) x $(dim(Hami)) with eig ...")
    tic()
    Hami.energies, Hami.eigenstates = eig(Hermitian(full(Hami.H)))
    toc()
    println()
    return Hami.energies, Hami.eigenstates
end

function eigs!(Hami::HaldaneSphereHamiLLL, nev::Int)
    @assert ishermitian(Hami.H)
    println("Diagonalizing Hami of size $(dim(Hami)) x $(dim(Hami)) with eigs, nev = $(nev) ...")
    tic()
    Hami.energies, Hami.eigenstates = eigs(Hami.H, which=:SR, nev=nev)
    toc()
    println()
    return Hami.energies, Hami.eigenstates
end

function organize_spectrum!(Hami::HaldaneSphereHamiLLL)
    println("Organizing spectrum according to quantum number L ...")
    tic()
    Hami.spectrum = zeros(Float64, (length(Hami.energies), 2))
    for j in 1:length(Hami.energies)
        Hami.spectrum[j, 1] = Lqn(Hami.eigenstates[:, j]' * Hami.setup.Lsquared * Hami.eigenstates[:, j])
        Hami.spectrum[j, 2] = Hami.energies[j]
    end
    toc()
    println()
    return Hami.spectrum
end

function entanglement_spectrum{T}(psi::Vector{T}, states::HaldaneSphereStatesListLLL, subsystemA::Vector{Int}, NA::Int, LzAvec::Vector{Float64})
    @assert length(psi) == dim(states)

    psi /= norm(psi) # normalize the wf

    IAlist = zeros(Int, length(psi))
    NAlist = zeros(Int, length(psi))
    LzAlist = zeros(Float64, length(psi))
    mzA = states.possible_mz[subsystemA]
    for i in 1:dim(states)
        nsA = states.nlist[i][subsystemA]
        IAlist[i] = I(nsA)
        NAlist[i] = sum(nsA)
        LzAlist[i] = vecdot(mzA, nsA)
    end

    println("Calculating entanglement spectrum for NA = $(NA), LzA = $(LzAvec) ...")
    tic()
    spec = zeros(Float64, (0, 2))
    for LzA in LzAvec

        # strip and sort
        strip_inds = intersect(find(NAlist .== NA), find(LzAlist .== LzA))
        IAlist_stripped = IAlist[strip_inds]
        sort_inds = sortperm(IAlist_stripped)
        IAlist_stripped = IAlist_stripped[sort_inds]
        # FIXME: check IB consistency
        psi_stripped = psi[strip_inds][sort_inds]

        # reshape psi into a matrix
        nc = length(find(IAlist_stripped .== IAlist_stripped[1]))
        nr = Int(length(psi_stripped) / nc)
        psi_matrix = reshape(psi_stripped, (nc, nr)).'

        # compute and diagonalize rhoA
        rhoA = psi_matrix * psi_matrix'
        evals = eigvals(rhoA)
        xi = -log.(evals[find(evals .> 0)])

        spec = vcat(spec, hcat(fill(LzA, length(xi)), sort(xi)))
    end
    toc()
    println()

    return spec
end

function entanglement_spectrum{T}(psi::Vector{T}, states::HaldaneSphereStatesListLLL, subsystemA::Vector{Int}, NA::Int, LzA::Float64)
    entanglement_spectrum(psi, states, subsystemA, NA, [LzA])
end

# matrix representation of c_alpha1p^dagger c_alpha2p^dagger c_alpha2 c_alpha1 in occupation number basis (very similiar to Hami building -- could probably reuse code..)
function four_point_operator(states::HaldaneSphereStatesListLLL, alpha1p::Int, alpha2p::Int, alpha1::Int, alpha2::Int)
    # make the following asserts to simplify the logic below (should be all we need with a little algebra)
    @assert alpha1p < alpha2p
    @assert alpha2 < alpha1
    @assert alpha1p + alpha2p == alpha1 + alpha2 # Lz conservation

    rows = Int[]
    cols = Int[]
    op_vals = Int[]
    @inbounds begin
    for j in 1:dim(states)
        # if j == 1 || mod(j, ceil(0.05dim(states))) == 0
        #     println("Working on c_$(alpha1p)^† c_$(alpha2p)^† c_$(alpha2) c_$(alpha1), column j = $(j): $(j)/$(dim(states)) = $(100j/dim(states))%")
        # end
        ns = states.nlist[j]
        if ns[alpha1] == 1 && ns[alpha2] == 1
            ns_new = copy(ns)
            ns_new[alpha1] = 0
            ns_new[alpha2] = 0
            if ns_new[alpha1p] == 0 && ns_new[alpha2p] == 0
                p1 = sum(ns_new[1:alpha1p-1])
                p2 = sum(ns_new[1:alpha2p-1])
                p3 = sum(ns[1:alpha1-1])
                p4 = sum(ns[1:alpha2-1])
                sgn = (-1)^(p1 + p2 + p3 + p4) # sign of matrix element from anticommutation relations
                ns_new[alpha1p] = 1
                ns_new[alpha2p] = 1
                i = searchsorted(states.Ilist, I(ns_new))[1] # binary search to figure out which row to populate
                push!(rows, i)
                push!(cols, j)
                push!(op_vals, sgn)
            end
        end
    end
    end
    op = sparse(rows, cols, op_vals, dim(states), dim(states))

    return op
end

function four_point_corrs{T}(psi::Vector{T}, states::HaldaneSphereStatesListLLL)
    @assert length(psi) == dim(states)

    Norb = get_Norb(states)

    println("Computing 4-point functions ...")
    tic()
    corr_mat = zeros(T, (Norb, Norb, Norb, Norb))
    for alpha2p in 1:Norb
        println("Working on alpha2p = $(alpha2p): $(alpha2p)/$(Norb)")
        for alpha1p in 1:alpha2p-1
            for alpha1 in 1:Norb
                for alpha2 in 1:alpha1-1
                    if alpha1p + alpha2p == alpha1 + alpha2 # Lz conservation [mz = alpha - (Q + 1) = alpha + const.]
                        corr_mat[alpha1p, alpha2p, alpha1, alpha2] = (psi' * four_point_operator(states, alpha1p, alpha2p, alpha1, alpha2) * psi) / (psi' * psi)
                    end
                end
            end
        end
    end
    toc()
    println()

    return corr_mat
end

function pair_correlation_function(psi::Vector, states::HaldaneSphereStatesListLLL,
                                   θ1::Vector{Float64}, ϕ1::Vector{Float64}, θ2::Vector{Float64}, ϕ2::Vector{Float64};
                                   n::Int = 0)
    @assert length(psi) == dim(states)
    @assert length(θ1) == length(ϕ1) == length(θ2) == length(ϕ2)
    @assert n >= 0

    ell = get_ell(states)
    Q = ell - n # to use orbitals corresponding to Landau level index n
    Norb = get_Norb(states)

    corr_mat = four_point_corrs(psi, states)

    println("Computing pair correlation function for $(length(θ1)) points ...")
    tic()
    g = zeros(Complex128, length(θ1))
    for i in 1:length(g)
        if i == 1 || mod(i, ceil(0.05length(g))) == 0
            println("Working on g[i], i = $(i): $(i)/$(length(g)) = $(100i/length(g))%")
        end
        for alpha2p in 1:Norb
            for alpha1p in 1:alpha2p-1
                for alpha1 in 1:Norb
                    for alpha2 in 1:alpha1-1
                        g1 = (conj(monopole_harmonic(Q, ell, m(alpha1p, ell), θ1[i], ϕ1[i]))
                            * conj(monopole_harmonic(Q, ell, m(alpha2p, ell), θ2[i], ϕ2[i]))
                            * monopole_harmonic(Q, ell, m(alpha2, ell), θ2[i], ϕ2[i])
                            * monopole_harmonic(Q, ell, m(alpha1, ell), θ1[i], ϕ1[i]))
                        g2 = (conj(monopole_harmonic(Q, ell, m(alpha1p, ell), θ1[i], ϕ1[i]))
                            * conj(monopole_harmonic(Q, ell, m(alpha2p, ell), θ2[i], ϕ2[i]))
                            * monopole_harmonic(Q, ell, m(alpha1, ell), θ2[i], ϕ2[i])
                            * monopole_harmonic(Q, ell, m(alpha2, ell), θ1[i], ϕ1[i]))
                        g3 = (conj(monopole_harmonic(Q, ell, m(alpha2p, ell), θ1[i], ϕ1[i]))
                            * conj(monopole_harmonic(Q, ell, m(alpha1p, ell), θ2[i], ϕ2[i]))
                            * monopole_harmonic(Q, ell, m(alpha2, ell), θ2[i], ϕ2[i])
                            * monopole_harmonic(Q, ell, m(alpha1, ell), θ1[i], ϕ1[i]))
                        g4 = (conj(monopole_harmonic(Q, ell, m(alpha2p, ell), θ1[i], ϕ1[i]))
                            * conj(monopole_harmonic(Q, ell, m(alpha1p, ell), θ2[i], ϕ2[i]))
                            * monopole_harmonic(Q, ell, m(alpha1, ell), θ2[i], ϕ2[i])
                            * monopole_harmonic(Q, ell, m(alpha2, ell), θ1[i], ϕ1[i]))
                        g[i] += (g1 - g2 - g3 + g4) * corr_mat[alpha1p, alpha2p, alpha1, alpha2]
                    end
                end
            end
        end
    end
    toc()
    println()

    return g
end

# evaulate along the equator using great circle (gc) distance (in units of ell0); g normalized to unity as r --> infty (divide by (density)^2)
# returns r/ell0 and g as column vectors
function pair_correlation_function_equator_gc(psi::Vector, states::HaldaneSphereStatesListLLL; pts::Int = 200, frac::Float64 = 0.5)
    ϕ = linspace(0, frac * 2pi, pts)
    g = pair_correlation_function(psi, states,
                                  fill(pi/2, length(ϕ)),
                                  collect(ϕ),
                                  fill(pi/2, length(ϕ)),
                                  fill(0., length(ϕ)),
                                  n=0)
    g /= (get_N(states)/(4pi))^2 # density = N/(4pi)
    println("Throwing out imaginary part of g: maximum(abs, imag(g)) = $(maximum(abs, imag(g)))")
    return hcat(sqrt(get_ell(states)) * ϕ, real(g))
end

max_stars_repo_name: mishmash/QHED

max_issues_repo_name: mishmash/QHED

max_forks_repo_name: mishmash/QHED


text:
using HampelFilter
using Test

@testset "HampelFilter.jl" begin
    # Write your tests here.
end