text:
using ProjAssistant
@quickactivate

# ----------------------------------------------------------------------------
## ARGS
using ArgParse

set = ArgParseSettings()
@add_arg_table! set begin
    "--new-dat"
        help = "ignore disk stored DAT"
        action = :store_true
    "--skip-me"
        help = "do not recompute MaxEnt part"
        action = :store_true
    "--skip-lp"
        help = "do not recompute LP part"
        action = :store_true
    "--dry-run"
        help = "do all but change nothing on disk"
        action = :store_true
end

if isinteractive()
    # Dev values
    new_dat = false
    skip_lp = false
    skip_me = false
    dry_run = false
else
    parsed_args = parse_args(set)
    new_dat = parsed_args["new-dat"]
    skip_lp = parsed_args["skip-lp"]
    skip_me = parsed_args["skip-me"]
    dry_run = parsed_args["dry-run"]
end

# ----------------------------------------------------------------------------
@time begin
    import Chemostat_EColi
    const ChE = Chemostat_EColi

    import Chemostat_Kayser2005
    const ChK = Chemostat_Kayser2005
    const KiJR = ChK.iJR904
    const Kd = ChK.KayserData

    import Chemostat_Nanchen2006
    const ChN = Chemostat_Nanchen2006
    const NiJR = ChN.iJR904
    const Nd = ChN.NanchenData

    import Chemostat_Folsom2014
    const ChF = Chemostat_Folsom2014
    const FiJR = ChF.iJR904
    const Fd = ChF.FolsomData

    import ChemostatPlots
    const ChP = ChemostatPlots

    import Chemostat
    import Chemostat.LP.MathProgBase
    const Ch = Chemostat
    const ChU = Ch.Utils
    const ChLP = Ch.LP
    using ProgressMeter
    using Base.Threads
    using SparseArrays

    using Serialization
    using MagicDicts

    using Statistics
end

## ----------------------------------------------------------------------------------
include("0.2.1_collect_commom.jl")
include("0.2.2_collect_LP.jl")
include("0.2.3_collect_ME.jl")

# ----------------------------------------------------------------------------------
let
    nths = nthreads()
    for (iJR, Ed, EXPS, EXCH_FLX_IDERS, INNER_FLX_IDERS) in [
        (NiJR, Nd, 1:17, ["GLC", "AC"], NiJR.load_inner_iders()),
        (KiJR, Kd, 5:12, ["GLC", "CO2", "O2", "AC", "NH4"], []),
        (FiJR, Fd, 1:4, ["GLC", "SUCC", "AC", "FORM"], []),
    ]
        src = string(nameof(Ed))

        # ---------------------------------------------------------------------------------
        DAT_FILE = ChE.DATfile(src)
        DAT = (!isfile(DAT_FILE) || new_dat) ? MagicDict() : ChE.load_DAT(src)

        DAT[:EXPS] = EXPS
        DAT[:EXCH_FLX_IDERS] = EXCH_FLX_IDERS
        DAT[:INNER_FLX_IDERS] = INNER_FLX_IDERS

        # LP methods
        DAT[:LP_METHODS] = [
            :FBA_Z_FIX_ADJ_MIN_VG_MIN_COST,
            :FBA_Z_FIX_ADJ_MIN_VG_MAX_COST,
            :FBA_Z_FIX_ADJ_MAX_VG_MIN_COST,
            :FBA_Z_FIX_ADJ_MAX_VG_MAX_COST,
            :FBA_Z_FIX_MIN_VG_MIN_COST,
            :FBA_Z_FIX_MIN_VG_MAX_COST,
            :FBA_Z_FIX_MAX_VG_MIN_COST,
            :FBA_Z_FIX_MAX_VG_MAX_COST,
            :FBA_MAX_Z_MIN_COST,
            :FBA_MAX_Z_MAX_COST,
            :FBA_Z_FIX_MIN_COST,
            :FBA_Z_FIX_MAX_COST,
            :FBA_Z_VG_FIX_MIN_COST,
            :FBA_Z_VG_FIX_MAX_COST,
            :FBA_Z_FIX_MAX_ATP_MIN_COST,
            :FBA_Z_FIX_MAX_ATP_MAX_COST
        ]

        # ME methods
        DAT[:ME_METHODS] = [
            :ME_MAX_POL,
            :ME_MAX_POL_B0,
            :ME_Z_EXPECTED_G_BOUNDED,
        ]

        DAT[:EXP_INDEPENDENT] = [:ME_MAX_POL_B0]

        # ---------------------------------------------------------------------------------
        collect_commom_DAT!(DAT, Ed, iJR)
        skip_lp || collect_lp_DAT!(DAT, Ed, iJR)
        skip_me || collect_ME_DAT!(DAT, iJR; nths)

        # ---------------------------------------------------------------------------------
        # saving
        dry_run || sdat(DAT, DAT_FILE)
    end

end;


text:
export save_ts_jld, save_ts_mat, parse_v
#=
methods to select which variables in workspace to save
=#

function parse_all(rhd)

    if rhd.save.save_full #parse voltage
        try
            v=parse_v(rhd.save.v)
        catch
            println("Couldn't parse voltage data!")
        end
    end
    if rhd.save.ttl_s
        try #TTL in
            ttlin=parse_ttl(rhd.save.ttl) #These are timestamps, polarity matters
        catch
            println("Could not parse TTL inputs")
        end

        #TTL out
        try
            ttlout=parse_ttl(rhd.save.ttl_out) #These are timestamps, polarity matters
        catch
            println("Could not parse TTL outputs")
        end
    end
    if rhd.save.adc_s
        try
            adcs=parse_adc(rhd.save.adc)
        catch
            println("Could not parse ADCs")
        end
    end
    if rhd.save.ts_s #Time stamps for sorting

    end
end

mutable struct Voltage_Header
    date_m::UInt8
    date_d::UInt8
    date_y::UInt16
    time_h::UInt8
    time_m::UInt8
    num_channels::UInt16
    samples_per_block::UInt16
end

Voltage_Header()=Voltage_Header(0,0,0,0,0,0,0)

if VERSION > v"0.7-"
    read(f,t::Type{UInt8},num) = Base.read(f,num)
    read(f,t::Type{UInt16},num) = reinterpret(UInt16,Base.read(f,num*2))
    read(f,t::Type{UInt32},num) = reinterpret(UInt32,Base.read(f,num*4))
    read(io::IO, T::Type) = Base.read(io::IO, T)
end

function read_v_header(fname="v.bin")

    f=open(fname,"r+")

    myheader=Voltage_Header()

    myheader.date_m=read(f,UInt8,1)[1]
    myheader.date_d=read(f,UInt8,1)[1]
    myheader.date_y=read(f,UInt16,1)[1]
    myheader.time_h=read(f,UInt8,1)[1]
    myheader.time_m=read(f,UInt8,1)[1]
    myheader.num_channels=read(f,UInt16,1)[1]
    myheader.samples_per_block=read(f,UInt16,1)[1]

    close(f)

    myheader
end

function prepare_lfp_header(rhd::RHD2000)

    if isfile(rhd.save.lfp)
        return
    end
    f=open(rhd.save.lfp,"a+")

    t=now()

    write(f,convert(UInt8,Dates.month(t)))
    write(f,convert(UInt8,Dates.day(t)))
    write(f,convert(UInt16,Dates.year(t)))
    write(f,convert(UInt8,Dates.hour(t)))
    write(f,convert(UInt8,Dates.minute(t)))
    write(f,convert(UInt16,size(rhd.v,2)))
    write(f,convert(UInt16,size(rhd.v,1)))

    close(f)

    nothing
end

function prepare_v_header(rhd::RHD2000)

    if isfile(rhd.save.v)
        return
    end

    f=open(rhd.save.v,"a+")

    t=now()

    write(f,convert(UInt8,Dates.month(t)))
    write(f,convert(UInt8,Dates.day(t)))
    write(f,convert(UInt16,Dates.year(t)))
    write(f,convert(UInt8,Dates.hour(t)))
    write(f,convert(UInt8,Dates.minute(t)))
    write(f,convert(UInt16,size(rhd.v,2)))
    write(f,convert(UInt16,size(rhd.v,1)))

    close(f)

    nothing
end

mutable struct Stamp_Header
    date_m::UInt8
    date_d::UInt8
    date_y::UInt16
    time_h::UInt8
    time_m::UInt8
    num_channels::UInt16
    sr::UInt32
end

Stamp_Header()=Stamp_Header(0,0,0,0,0,0,0)

function read_stamp_header(fname="ts.bin")

    f=open(fname,"r+")

    myheader=Stamp_Header()

    myheader.date_m=read(f,UInt8,1)[1]
    myheader.date_d=read(f,UInt8,1)[1]
    myheader.date_y=read(f,UInt16,1)[1]
    myheader.time_h=read(f,UInt8,1)[1]
    myheader.time_m=read(f,UInt8,1)[1]
    myheader.num_channels=read(f,UInt16,1)[1]
    myheader.sr=read(f,UInt32,1)[1]

    close(f)

    myheader
end

function prepare_stamp_header(rhd::RHD2000)

    if isfile(rhd.save.ts)
        return
    end

    f=open(rhd.save.ts,"a+")

    t=now()

    write(f,convert(UInt8,Dates.month(t)))
    write(f,convert(UInt8,Dates.day(t)))
    write(f,convert(UInt16,Dates.year(t)))
    write(f,convert(UInt8,Dates.hour(t)))
    write(f,convert(UInt8,Dates.minute(t)))
    write(f,convert(UInt16,size(rhd.v,2)))
    write(f,convert(UInt32,rhd.sr))

    close(f)

    nothing
end

mutable struct ADC_Header
    date_m::UInt8
    date_d::UInt8
    date_y::UInt16
    time_h::UInt8
    time_m::UInt8
    num_channels::UInt16
    samples_per_block::UInt16
end

ADC_Header()=ADC_Header(0,0,0,0,0,0,0)

function read_adc_header(fname="adc.bin")

    f=open(fname,"r+")

    myheader=ADC_Header()

    myheader.date_m=read(f,UInt8,1)[1]
    myheader.date_d=read(f,UInt8,1)[1]
    myheader.date_y=read(f,UInt16,1)[1]
    myheader.time_h=read(f,UInt8,1)[1]
    myheader.time_m=read(f,UInt8,1)[1]
    myheader.num_channels=read(f,UInt16,1)[1]
    myheader.samples_per_block=read(f,UInt16,1)[1]

    close(f)

    myheader
end

function prepare_adc_header(rhd::RHD2000)

    if isfile(rhd.save.adc)
        return
    end

    f=open(rhd.save.adc,"a+")

    t=now()

    write(f,convert(UInt8,Dates.month(t)))
    write(f,convert(UInt8,Dates.day(t)))
    write(f,convert(UInt16,Dates.year(t)))
    write(f,convert(UInt8,Dates.hour(t)))
    write(f,convert(UInt8,Dates.minute(t)))
    write(f,convert(UInt16,8))
    write(f,convert(UInt16,size(rhd.v,1)))

    close(f)

    nothing
end

mutable struct TTL_Header
    date_m::UInt8
    date_d::UInt8
    date_y::UInt16
    time_h::UInt8
    time_m::UInt8
    num_channels::UInt16
    samples_per_block::UInt16
    sr::UInt32
end

TTL_Header()=TTL_Header(0,0,0,0,0,0,0,0)

function read_ttl_header(fname="ttl.bin")

    f=open(fname,"r+")

    myheader=TTL_Header()

    myheader.date_m=read(f,UInt8,1)[1]
    myheader.date_d=read(f,UInt8,1)[1]
    myheader.date_y=read(f,UInt16,1)[1]
    myheader.time_h=read(f,UInt8,1)[1]
    myheader.time_m=read(f,UInt8,1)[1]
    myheader.num_channels=read(f,UInt16,1)[1]
    myheader.samples_per_block=read(f,UInt16,1)[1]
    myheader.sr=read(f,UInt32,1)[1]

    close(f)

    myheader
end

function prepare_ttl_header(rhd::RHD2000,savepath=rhd.save.ttl)

    if isfile(savepath)
        return
    end

    f=open(savepath,"a+")

    t=now()

    write(f,convert(UInt8,Dates.month(t)))
    write(f,convert(UInt8,Dates.day(t)))
    write(f,convert(UInt16,Dates.year(t)))
    write(f,convert(UInt8,Dates.hour(t)))
    write(f,convert(UInt8,Dates.minute(t)))
    write(f,convert(UInt16,16))
    write(f,convert(UInt16,size(rhd.v,1)))
    write(f,convert(UInt32,rhd.sr))

    close(f)

    nothing
end

#=
Voltage
=#

function parse_v(fname="v.bin")

    myheader=read_v_header(fname)

    f=open(fname, "r+")

    seekend(f)
    l=position(f)-10
    v=zeros(Int16,div(l,2*myheader.num_channels),myheader.num_channels)

    seek(f,10)

    count=0
    while eof(f)==false
        for i=1:myheader.num_channels
            for j=1:myheader.samples_per_block
                v[count+j,i]=read(f,Int16)
            end
        end
        count+=myheader.samples_per_block
    end
    close(f)

    v
end

read_single_v(rhd::RHD2000,chan) = read_single_v(rhd.save.v,chan)

function read_single_v(fname::String,chan)

    myheader=read_v_header(fname)

    f=open(fname,"r+")

    seekend(f)

    l=position(f)-10
    v=zeros(Int16,div(l,2*myheader.num_channels))

    seek(f,10)

    count=0
    while eof(f)==false
        for i=1:myheader.num_channels
            if i==chan
                for j=1:myheader.samples_per_block
                    v[count+j]=read(f,Int16)
                end
            else
                for j=1:myheader.samples_per_block
                    read(f,Int16)
                end
            end
        end
        count+=myheader.samples_per_block
    end
    close(f)

    v
end

function save_v_mat(in_name="v.bin",out_name="v.mat")

    v=parse_v(in_name)

    file = matopen(out_name, "w")
    write(file, "v", v)
    close(file)

    nothing
end

function save_v_jld(in_name="v.bin",out_name="v.jld")

    v=parse_v(in_name)

    file = jldopen(out_name, "w")
    write(file, "v", v)
    close(file)

    nothing
end

#=
Time Stamps
=#

function parse_ts(fname="ts.bin")

    myheader=read_stamp_header(fname)

    ss=[Array(Spike,0) for i=1:myheader.num_channels]

    numcells=zeros(Int64,myheader.num_channels)

    f=open(fname, "r+")

    seek(f,12)

    while eof(f)==false

        t=read(f,UInt32,1)[1]

        for j=1:myheader.num_channels
            chan=read(f,UInt16,1)[1] #channel
            num=read(f,UInt16,1)[1] #Number of upcoming spikes
            for i=1:num
                myss=read(f,Int64,2) #time stamps
                clus=read(f,UInt8,1)[1] #cluster
                if clus>numcells[j]
                    numcells[j]=clus
                end
                push!(ss[j],Spike((t+myss[1]:t+myss[2]),clus))
            end
        end

    end

    close(f)

    (ss,numcells)
end

function get_ts_dict(ss::Array{Array{Spike,1},1},numcells::Array{Int64,1},sr=30000,tmin=0.0)

    spikes=Dict{String,Array{Float64,1}}()

    for i=1:length(numcells)
        for j=1:numcells[i]
            myname=string("s",i,"_",j)

            myspikes=zeros(Float64,0)

            for k=1:length(ss[i])
                if (ss[i][k].inds.start/sr>tmin) &&(ss[i][k].id==j)
                    push!(myspikes,ss[i][k].inds.start/sr)
                end
            end
            spikes[myname]=myspikes
        end
    end
    spikes
end

function save_ts_mat(in_name="ts.bin",out_name="ts.mat")

    myheader=read_stamp_header(in_name)

    (ss,numcells)=parse_ts(in_name)

    spikes=get_ts_dict(ss,numcells,myheader.sr)

    file = matopen(out_name, "w")
    write(file, "spikes", spikes)
    close(file)

    spikes
end

function save_ts_jld(in_name="ts.bin",out_name="ts.jld")

    myheader=read_stamp_header(in_name)

    (ss,numcells)=parse_ts(in_name)

    spikes=get_ts_dict(ss,numcells,myheader.sr)

    file = jldopen(out_name, "w")
    write(file, "spikes", spikes)
    close(file)

    spikes
end

#=
ADC Signals
=#

function parse_adc(fname="adc.bin")

    myheader=read_adc_header(fname)

    f=open(fname, "r+")

    seekend(f)
    l=position(f)-10
    adc=zeros(UInt16,div(l,2*myheader.num_channels),myheader.num_channels)

    seek(f,10)

    count=0
    while eof(f)==false
        for i=1:myheader.num_channels
            for j=1:myheader.samples_per_block
                adc[count+j,i]=read(f,UInt16)
            end
        end
        count+=myheader.samples_per_block
    end
    close(f)

    adc
end

#=
Events
=#

function parse_ttl(fname="ttl.bin",polarity=trues(16))

    myheader=read_ttl_header(fname)

    positive_event=trues(16);

    for i=1:length(polarity)
        positive_event[i]=polarity[i]
    end

    ttl_times=[zeros(Int64,0) for i=1:myheader.num_channels]

    f=open(fname, "r+")

    seek(f,14) #start at second time step

    x_p=read(f,UInt16)

    count=2
    while eof(f)==false

        x=read(f,UInt16)

        for i=1:myheader.num_channels
            y=x&(2^(i-1))
            if (y>0)==(positive_event[i])
                #y_p=x_p&(2^(i-1))
                #if y_p==0
                    push!(ttl_times[i],count)
                #end
            end
        end
        #x_p=x
        count+=1
    end
    close(f)

    ttl_times
end

#=
EXPORT TO PLEXON

Methods to convert data waveforms to PLX so it can be used with offline sorter.
Thanks to Simon Kornblith for linking to description of PLX data structures here:
http://hardcarve.com/wikipic/PlexonDataFileStructureDocumentation.pdf
=#

mutable struct PL_FileHeader
	MagicNumber::UInt32
	Version::Int32
    	Comment::Array{UInt8,1}
	ADFrequency::Int32
	NumDSPChannels::Int32
	NumEventChannels::Int32
	NumSlowChannels::Int32
	NumPointsWave::Int32
	NumPointsPreThr::Int32
	Year::Int32
	Month::Int32
	Day::Int32
	Hour::Int32
	Minute::Int32
	Second::Int32
	FastRead::Int32
	WaveformFreq::Int32
	LastTimestamp::Float64
	Trodalness::UInt8
	DataTrodalness::UInt8
	BitsPerSpikeSample::UInt8
	BitsPerSlowSample::UInt8
	SpikeMaxMagnitudeMV::UInt16
	SlowMaxMagnitudeMV::UInt16
	SpikePreAmpGain::UInt16
    	Padding::Array{UInt8,1}
    	TSCounts::Array{Int32,2}
    	WFCounts::Array{Int32,2}
    	EVCounts::Array{Int32,1}
end

function PL_FileHeader(sr,num_chan,num_point,pre_t,t_end,tscounts,event_chan,evcounts)
    mn=0x58454c50
    ver=105
    com=zeros(UInt8,128)
    pad=zeros(UInt8,46)

    PL_FileHeader(mn,ver,com,sr,num_chan,event_chan,0,num_point,pre_t,1,1,1,1,1,1,0,0,t_end,0x01,0x01,16,12,3000,5000,0x03e8,pad,tscounts,tscounts,evcounts)
end

mutable struct PL_ChanHeader
    	Name::Array{UInt8,1}
    	SIGName::Array{UInt8,1}
	Channel::Int32
	WFRate::Int32
	SIG::Int32
	Ref::Int32
	Gain::Int32
	Filter::Int32
	Threshold::Int32
	Method::Int32
	NUnits::Int32
    	Template::Array{Int16,2}
    	Fit::Array{Int32,1}
	SortWidth::Int32
    	Boxes::Array{Int16,3}
	SortBeg::Int32
    	Comment::Array{UInt8,1}
    	Padding::Array{Int32,1}
end

function PL_ChanHeader(num,units)

    if num<10
        c=string("sig00",num)
    elseif num<100
        c=string("sig0",num)
    else
        c=string("sig",num)
    end
    myname=zeros(UInt8,32)
    for i=1:length(c)
        myname[i]=convert(UInt8,c[i])
    end
    templates=zeros(Int16,5,64)
    fits=zeros(Int32,5)
    boxes=zeros(Int16,5,2,4)
    com=zeros(UInt8,128)
    pad=zeros(Int32,11)

    PL_ChanHeader(myname,myname,num,10,num,0,16,0,0,1,units,templates,fits,0,boxes,0,com,pad)
end

mutable struct PL_EventHeader
    Name::Array{UInt8,1}
    Channel::Int32
    Comment::Array{UInt8,1}
    Padding::Array{Int32,1}
end

function PL_EventHeader(num)
    if num<10
        c=string("event0",num)
    else
        c=string("event",num)
    end
    myname=zeros(UInt8,32)
    for i=1:length(c)
        myname[i]=convert(UInt8,c[i])
    end
    PL_EventHeader(myname,num,zeros(UInt8,128),zeros(Int32,33))
end

mutable struct PL_DataBlockHeader
    Type::Int16
    UpperByteOf5ByteTimestamp::UInt16
    TimeStamp::UInt32
    Channel::Int16
    Unit::Int16
    NumberOfWaveForms::Int16
    NumberOfWordsInWaveform::Int16
end

function PL_DataBlockHeader(mytype,t,num,unit,wave_size)
    if mytype==1
        PL_DataBlockHeader(1,0,t,num,unit,1,wave_size)
    elseif mytype==4
        PL_DataBlockHeader(4,0,t,num,0,0,0)
    end
end

function write_plex(out_name::AbstractString,vname="v.bin",tsname="ts.bin"; ttl_parse=false,ttlname="ttl.bin",tmin=0)

    v_header=read_v_header(vname)
    ts_header=read_stamp_header(tsname)

    sample_size=v_header.samples_per_block
    sr=ts_header.sr

    (ss,numcells)=parse_ts(tsname)

    num_channel=length(ss)

    spikes=get_ts_dict(ss,numcells,tmin,sr)

    tscounts=zeros(Int32,130,5)

	#For plexon, 0 is unsorted, 1-4 is sorted
	#Intan will save starting at 1 for unsorted, and all numbers after as units
	#This array only logs sorted units. 0 entry (1 in Julia) is unused
    for i=1:length(ss)
        for j=1:length(ss[i])
	if ss[i][j].id>1
            tscounts[i+1,ss[i][j].id]+=1
	end
	end
    end

    samples_per_wave=convert(Int16,length(ss[1][1].inds))

	if ttl_parse
    		myttl=parse_ttl(ttlname)
	end
    evcounts=zeros(Int32,512)

	if ttl_parse
    for i=1:16
        evcounts[i]=length(myttl[i])
    end
	end

    f_out=open(out_name,"a+")

    #First is File Header
    file_header=PL_FileHeader(sr,length(ss),samples_per_wave,24,ss[1][end].inds.stop/sr,tscounts,16,evcounts)

    for i=1:length(fieldnames(file_header))
        write(f_out,getfield(file_header,i))
    end

    #Then Spike Channel Headers
    for i=1:length(ss)
        chan_header=PL_ChanHeader(i,numcells[i])
        for j=1:length(fieldnames(chan_header))
            write(f_out,getfield(chan_header,j))
        end
    end

    #Finally is event channel headers
    for i=1:16
        event_header=PL_EventHeader(i)
        for j=1:length(fieldnames(event_header))
            write(f_out,getfield(event_header,j))
        end
    end

    #This will suck for big files
    #should read in one channel voltage at a time

    #Write data blocks for spike waveforms
    myv=zeros(Int16,samples_per_wave)
    for i=1:length(ss)
	v=read_single_v(vname,i)
        for j=1:length(ss[i])
            header=PL_DataBlockHeader(1,ss[i][j].inds.start,i,ss[i][j].id-1,samples_per_wave)
            myind=ss[i][j].inds.start

            if myind+samples_per_wave-1 < size(v,1)
                count=1
                for k=myind:(myind+samples_per_wave-1)
                    myv[count]=v[k]
                    count+=1
                end
                for k=1:length(fieldnames(header))
                    write(f_out,getfield(header,k))
                end
                write(f_out,myv)
            end
        end
    end

    #Write Event blocks
	if ttl_parse
    for i=1:16
        for j=1:length(myttl[i])
            header=PL_DataBlockHeader(4,myttl[i][j],i,0,0)
            for k=1:length(fieldnames(header))
                write(f_out,getfield(header,k))
            end
        end
    end
	end

    close(f_out)

    nothing
end

function save_config_cb(widget::Ptr,user_data::Tuple{Gui_Handles,R,Array{S,1}}) where {R<:RHD2000,S<:Sorting}

    han, rhd, s = user_data

    filepath=save_dialog("Save configuration",han.win)

    if filepath != ""

        if filepath[end-3:end]==".jld"
        else
            filepath=string(filepath,".jld")
        end

        file = jldopen(filepath, "w")

        write(file, "Gain", han.scale)
        write(file, "Offset", han.offset)
        write(file, "Sorting", s)
        write(file, "total_clus",han.total_clus)
        write(file, "Enabled", han.enabled)
        write(file, "Reference",rhd.refs)

        close(file)

    end

    nothing
end

function export_plex_cb(widget::Ptr,user_data::Tuple{Gui_Handles,RHD2000})

    han, rhd = user_data

    write_plex(save_dialog("Export to Plex",han.win),rhd.save.v,rhd.save.ts,rhd.save.ttl)

    nothing
end

function export_jld_cb(widget::Ptr,user_data::Tuple{Gui_Handles,RHD2000})

    han,rhd = user_data

    #get prefix for exported data name

    #find out what to export

    #call appropriate save functions based on above

    nothing
end

function export_mat_cb(widget::Ptr,user_data::Tuple{Gui_Handles,RHD2000})

    han,rhd = user_data

    nothing
end

function load_config_cb(widget::Ptr,user_data::Tuple{Gui_Handles,R,Array{S,1}}) where {R<:RHD2000,S<:Sorting}

    han, rhd, s = user_data

    filepath=open_dialog("Load Configuration",han.win)

    if filepath != ""

        c = jldopen(filepath, "r") do file
            g=read(file,"Gain")

            for i=1:length(g)
                han.scale[i]=g[i]
            end

            o=read(file,"Offset")

            for i=1:length(o)
                han.offset[i]=o[i]
            end

            s_saved=read(file,"Sorting")

            for i=1:length(s)
                s[i]=s_saved[i]
            end

            total_clus=read(file,"total_clus")
            for i=1:length(total_clus)
                han.total_clus[i]=total_clus[i]
            end

            e=read(file,"Enabled")

            for i=1:length(e)
                han.enabled[i]=e[i]
            end

            refs=read(file,"Reference")

            for i=1:length(refs)
                rhd.refs[i]=refs[i]
            end
        end

        update_treeview(han.sc)

        update_ref(rhd,han)

    end

    nothing
end

function load_backup_cb(widget::Ptr,user_data::Tuple{Gui_Handles,R,Array{S,1}}) where {R<:RHD2000,S<:Sorting}

    han, rhd, s = user_data

    filepath=open_dialog("Select Backup File",han.win)

    if filepath != ""

        filepath=filepath[1:(end-10)]
        println(filepath)
        #Enabled

        f=open(string(filepath,"enabled.bin"),"r")

        seekend(f)
        chan_num=(position(f)<size(rhd.v,2)) ? position(f) : size(rhd.v,2)

        seekstart(f)
        for i=1:(chan_num)
            han.enabled[i]=read(f,UInt8)
        end

        close(f)

        #Reference
        if isfile(string(filepath,"ref.bin"))
            f=open(string(filepath,"ref.bin"),"r")
            for i=1:(chan_num)
                rhd.refs[i]=read(f,Int64)
            end
            close(f)
        end

        #Gain
        for i=1:chan_num
            if isfile(string(filepath,"gain/",i,".bin"))
                f=open(string(filepath,"gain/",i,".bin"),"r")
                han.scale[i,1]=read(f,Float64)
                han.scale[i,2]=han.scale[i,1]*.2
                close(f)
            end
        end

        #Thres
        for i=1:chan_num
            if isfile(string(filepath,"thres/",i,".bin"))
                f=open(string(filepath,"thres/",i,".bin"),"r")
                s[i].thres=read(f,Float64)
                close(f)
            end
        end

        for i=1:chan_num
            if isfile(string(filepath,"cluster/",i,".bin"))
                f=open(string(filepath,"cluster/",i,".bin"),"r")
                for j=1:length(han.sc.temp.templates)
                    s[i].c.templates[j]=read(f,Float64)
                end
                for j=1:length(han.sc.temp.sig_min)
                    s[i].c.sig_min[j]=read(f,Float64)
                end
                for j=1:length(han.sc.temp.sig_max)
                    s[i].c.sig_max[j]=read(f,Float64)
                end
                s[i].c.misses=read(f,Int64)
                s[i].c.num=read(f,Int64)
                for j=1:length(han.sc.temp.tol)
                    s[i].c.tol[j]=read(f,Float64)
                end
                han.total_clus[i]=s[i].c.num
                close(f)
            end
        end

        get_thres(han,s)
        get_cluster(han,s)
        SpikeSorting.update_treeview(han.sc)
        SpikeSorting.select_unit(han.sc)
        SpikeSorting.draw_templates(han.sc)
        update_ref(rhd,han)
    end
    nothing
end

ext: jl


text:
function cuITensor(::Type{T},inds::IndexSet) where {T<:Number}
    return ITensor(Dense{float(T)}(CuArrays.zeros(float(T),dim(inds))), inds)
end
cuITensor(::Type{T},inds::Index...) where {T<:Number} = cuITensor(T,IndexSet(inds...))

cuITensor(is::IndexSet)   = cuITensor(Float64,is)
cuITensor(inds::Index...) = cuITensor(IndexSet(inds...))

cuITensor() = ITensor()
function cuITensor(x::S, inds::IndexSet{N}) where {S<:Number, N}
    dat = CuVector{float(S)}(undef, dim(inds))
    fill!(dat, float(x))
    ITensor(Dense{S}(dat), inds)
end
cuITensor(x::S, inds::Index...) where {S<:Number} = cuITensor(x,IndexSet(inds...))

#TODO: check that the size of the Array matches the Index dimensions
function cuITensor(A::Array{S},inds::IndexSet) where {S<:Number}
    return ITensor(Dense(CuArray{S}(A)), inds)
end
function cuITensor(A::CuArray{S},inds::IndexSet) where {S<:Number}
    return ITensor(Dense(A), inds)
end
cuITensor(A::Array{S},   inds::Index...) where {S<:Number} = cuITensor(A,IndexSet(inds...))
cuITensor(A::CuArray{S}, inds::Index...) where {S<:Number} = cuITensor(A,IndexSet(inds...))
cuITensor(A::ITensor) = store(tensor(A)) isa ITensors.Empty ? cuITensor(A.inds) : cuITensor(data(store(tensor(A))), A.inds)

function Base.collect(A::ITensor)
    typeof(A.store.data) <: CuArray && return ITensor(collect(A.store), A.inds)
    return A
end

function randomCuITensor(::Type{S},inds::IndexSet) where {S<:Real}
  T = cuITensor(S,inds)
  randn!(T)
  return T
end
function randomCuITensor(::Type{S},inds::IndexSet) where {S<:Complex}
  Tr = cuITensor(real(S),inds)
  randn!(Tr)
  Ti = cuITensor(real(S),inds)
  randn!(Ti)
  return complex(Tr) + im.*complex(Ti)
end
randomCuITensor(::Type{S},inds::Index...) where {S<:Number} = randomCuITensor(S,IndexSet(inds...))
randomCuITensor(inds::IndexSet) = randomCuITensor(Float64,inds)
randomCuITensor(inds::Index...) = randomCuITensor(Float64,IndexSet(inds...))

CuArrays.CuArray(T::ITensor) = CuArray(tensor(T))

CuArrays.CuArray(T::ITensor,ninds::Index...) = storage_convert(CuArray,store(T),inds(T),IndexSet(ninds))

CuArrays.CuMatrix(A::ITensor) = CuArray(A)

function CuVector(A::ITensor)
  if ndims(A) != 1
    throw(DimensionMismatch("Vector() expected a 1-index ITensor"))
  end
  return CuArray(A)
end

function CuArray(T::ITensor{N},is::Vararg{Index,N}) where {N}
  perm = getperm(inds(T),is)
  return CuArray(permutedims(tensor(T),perm))
end

function CuMatrix(T::ITensor{N},i1::Index,i2::Index) where {N}
  N≠2 && throw(DimensionMismatch("ITensor must be order 2 to convert to a Matrix"))
  return CuArray(T,i1,i2)
end


text:
"""
Construct a vector of repeated values, just like R's `rep` function.
We do not have a `length_out` argument at this time.

# Examples

```jldoctest
julia> rep(["Go", "Fight", "Win"], times=2)
6-element Array{String,1}:
 "Go"
 "Fight"
 "Win"
 "Go"
 "Fight"
 "Win"
```

```jldoctest
julia> rep(["A", "B", "C"], each=3)
9-element Array{String,1}:
 "A"
 "A"
 "A"
 "B"
 "B"
 "B"
 "C"
 "C"
 "C"
```
"""
function rep(x::Union{Any,Vector}; each::Union{Int,Vector{Int}} = ones(Int,length(x)), times::Int = 1)
  if !isa(x,Vector)
    x = [ x ]
  end
  if isa(each,Int)
    each = [ each for i in eachindex(x) ]
  end
  length(x) != length(each) && throw(ArgumentError("If the arguemnt 'each' is not a scalar, it must have the same length as 'x'."))
  length_out = sum(each * times)
  rval = similar(x,length_out)
  index = 0
  for i in 1:times
    for j in eachindex(x)
      for k in 1:each[j]
        index += 1
        rval[index] = x[j]
      end
    end
  end
  return(rval)
end

"""
The four argument version substitutes customized ordering for a hard-coded '<'.
"""
function Base.searchsortedfirst(v::AbstractVector, x, lo::Int, hi::Int)
    lo = lo - 1
    hi = min(length(v), hi) + 1
    @inbounds while lo < hi-1
        m = (lo+hi)>>>1
        if v[m] < x
            lo = m
        else
            hi = m
        end
    end
    return hi
end

"""
The method for two vectors is like R's findinterval.
"""
Base.searchsortedfirst(v::AbstractVector, x::AbstractVector) = searchsortedfirst(v, x, 1, length(v))

function searchsortedfirst(v::AbstractVector, x::AbstractVector, lo::Int, hi::Int)
    indices = similar(x)
    min = lo - 1
    max = hi + 1
    @inbounds for (i,query) in enumerate(x)
        hi = hi + 1 # 2X speedup with this *inside* the loop for sorted x
        # unsorted x, restart left side
        if lo <= min || query <= v[lo]
            lo = min
        end
        if hi >= max || query >= v[hi]
            hi = max
        end
        # binary search for the exact bin
        while hi - lo > 1
            m = (lo+hi)>>>1
            if query > v[m]
                lo = m
            else
                hi = m
            end
        end
        indices[i] = hi
    end
    return(indices)
end


text:
@enum GlfwClientApi begin
    GlfwClientApi_Unknown
    GlfwClientApi_OpenGL
    GlfwClientApi_Vulkan
end

const GLFW_HAS_WINDOW_TOPMOST = true
const GLFW_HAS_WINDOW_HOVERED = true
const GLFW_HAS_WINDOW_ALPHA = true
const GLFW_HAS_PER_MONITOR_DPI = true
const GLFW_HAS_VULKAN = true
const GLFW_HAS_FOCUS_WINDOW = true
const GLFW_HAS_FOCUS_ON_SHOW = true
const GLFW_HAS_MONITOR_WORK_AREA = true
const GLFW_HAS_OSX_WINDOW_POS_FIX = true
const GLFW_HAS_NEW_CURSORS = false          # pending GLFW_jll@3.4
const GLFW_HAS_MOUSE_PASSTHROUGH = false    # pending GLFW_jll@3.4

const GLFW_WINDOW_DEFAULT_SIZE = Ref{Tuple{Int,Int}}((640,480))
const GLFW_WINDOW_DEFAULT_TITLE = Ref{String}("IMGUI")

const GLFW_BACKEND_PLATFORM_NAME = "julia_imgui_glfw"


text:
"""
Simulation of airspace consisting of `n` aircraft equipped with collision avoidance systems.
"""
module SimpleACAS

using AdaStress
using Dates
using Distributions
using IterTools
using LinearAlgebra
using Plots
using ProgressMeter
using Random

export visualize

include("utils.jl")
include("aircraft.jl")
include("encounter.jl")
include("cas.jl")
include("response.jl")
include("metrics.jl")
include("simulator.jl")
include("logging.jl")
include("visualization.jl")

end


text:
# __BEGIN_LICENSE__
#
# ThreeDeconv.jl
#
# Copyright (c) 2018, Stanford University
#
# All rights reserved.
#
# Redistribution and use in source and binary forms for academic and other
# non-commercial purposes with or without modification, are permitted provided
# that the following conditions are met:
#
# * Redistributions of source code, including modified source code, must retain
#   the above copyright notice, this list of conditions and the following
#   disclaimer.
#
# * Redistributions in binary form or a modified form of the source code must
#   reproduce the above copyright notice, this list of conditions and the
#   following disclaimer in the documentation and/or other materials provided with
#   the distribution.
#
# * Neither the name of The Leland Stanford Junior University, any of its
#   trademarks, the names of its employees, nor contributors to the source code
#   may be used to endorse or promote products derived from this software without
#   specific prior written permission.
#
# * Where a modified version of the source code is redistributed publicly in
#   source or binary forms, the modified source code must be published in a freely
#   accessible manner, or otherwise redistributed at no charge to anyone
#   requesting a copy of the modified source code, subject to the same terms as
#   this agreement.
#
# THIS SOFTWARE IS PROVIDED BY THE TRUSTEES OF THE LELAND STANFORD JUNIOR
# UNIVERSITY "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE LELAND STANFORD JUNIOR
# UNIVERSITY OR ITS TRUSTEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
# BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
# IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
# OF SUCH DAMAGE.
#
# __END_LICENSE__

using FFTW, Statistics
import Base.copy!, LinearAlgebra.norm

function circconv(
    x::AbstractMatrix{T}, h1::AbstractVector{T}, h2::AbstractVector{T})::AbstractMatrix{T} where T<:Real
    N1, N2 = size(x)
    X = rfft(x)
    if h1 == [T(1)]
        h_circ2 = padcircshift(h2, N2)
        H_circ2 = fft(h_circ2)
        out = irfft( X  .* transpose(H_circ2), size(x,1))
    elseif h2 == [T(1)]
        h_circ1 = padcircshift(h1, N1)
        H_circ1 = rfft(h_circ1)
        out = irfft( X .* H_circ1, size(x,1))
    else
        h_circ1 = padcircshift(h1, N1)
        H_circ1 = rfft(h_circ1)
        h_circ2 = padcircshift(h2, N2)
        H_circ2 = fft(h_circ2)
        out = irfft((X .* H_circ1) .* transpose(H_circ2), size(x,1))
    end
    return out
end

function circconv(
    x::AbstractMatrix{T}, h::AbstractMatrix{T}) where T<:Real
    N1, N2 = size(x)
    k1, k2 = size(h)
    p1 = mod(k1,2) == 0 ? div(k1,2) : div(k1-1,2)
    p2 = mod(k2,2) == 0 ? div(k2,2) : div(k2-1,2)
    h_pad = zeros(T, N1, N2)
    h_pad[1:k1, 1:k2] = h
    h_circ = circshift(h_pad, (-p1, -p2))
    return irfft(rfft(x) .* rfft(h_circ), size(x,1))
end

function circconv(x::AbstractMatrix{T}, h::AbstractVector{T}) where T<:Real
    N1, N2 = size(x)
    X = rfft(x)
    h_circ1 = padcircshift(h, N1)
    H_circ1 = rfft(h_circ1)
    h_circ2 = padcircshift(h, N2)
    H_circ2 = fft(h_circ2)
    return irfft((X .* H_circ1) .* transpose(H_circ2), size(x,1))
end

circconv1dx(x::Matrix{T}, h::Vector{T}) where T = circconv(x, [1.0], h)
circconv1dy(x::Matrix{T}, h::Vector{T}) where T = circconv(x, h, [1.0])

function padcircshift(h::AbstractVector{T}, N::Integer) where T<:Real
  k = length(h)
  @assert N > k
  h_pad = zeros(T, N)
  h_pad[1:k] = h
  p = mod(k,2) == 0 ? div(k,2) : div(k-1,2)
  return circshift(h_pad, -p)
end

mad(v::Vector{T}, κ::Float64 = 0.6745) where T= 1.0 / κ * median(abs.(v))

function madBiasFactor(n::Int)::Float64
    lut = [0.798, 0.732, 0.712, 0.702, 0.696, 0.693, 0.690, 0.688, 0.686, 0.685]
    if n <= 20
        out = lut[div(n-1,2) + 1]
    else
        out = 1 / (5.0 * n) + 0.6745
    end
    return out
end

function center_zeropad2!(
    padded_u::AbstractArray{T,2}, u::AbstractArray{T,2}) where T
    M, N = size(u)
    m1 = isodd(M) ? div(M+1,2) : div(M,2)
    m2 = M - m1
    n1 = isodd(N) ? div(N+1,2) : div(N,2)
    n2 = N - n1

    padded_u[1:m1, 1:n1] = u[1:m1, 1:n1]
    padded_u[1:m1, end-n2+1:end] = u[1:m1, end-n2+1:end]
    padded_u[end-m2+1:end, 1:n1] = u[end-m2+1:end, 1:n1]
    padded_u[end-m2+1:end, end-n2+1:end] = u[end-m2+1:end, end-n2+1:end]
end

function center_zeropad2(u::AbstractArray{T,2}, padded_shape::NTuple{2,Int}) where T

    padded_u = zeros(eltype(u), padded_shape)
    center_zeropad2!(padded_u, u)

    return padded_u
end

function center_zeropad2(
    u::AbstractArray{T,3}, padded_shape::NTuple{3,Int}) where T
    @assert size(u,3) == padded_shape[3]
    K = size(u,3)
    padded_u = zeros(eltype(u), padded_shape)

    for k in 1:K
        center_zeropad2!(view(padded_u,:,:,k), view(u,:,:,k))
    end

    return padded_u
end

function center_zeropad3!(
    padded_u::AbstractArray{T,3}, u::AbstractArray{T,3}) where T
    I, J, K = size(u)
    i1 = isodd(I) ? div(I+1,2) : div(I,2)
    i2 = I - i1
    j1 = isodd(J) ? div(J+1,2) : div(J,2)
    j2 = J - j1
    k1 = isodd(K) ? div(K+1,2) : div(K,2)
    k2 = K - k1

    padded_u[1:i1, 1:j1, 1:k1] = u[1:i1, 1:j1, 1:k1]
    padded_u[1:i1, end-j2+1:end, 1:k1] = u[1:i1, end-j2+1:end, 1:k1]
    padded_u[end-i2+1:end, 1:j1, 1:k1] = u[end-i2+1:end, 1:j1, 1:k1]
    padded_u[end-i2+1:end, end-j2+1:end, 1:k1] = u[end-i2+1:end, end-j2+1:end, 1:k1]
    padded_u[1:i1, 1:j1, end-k2+1:end] = u[1:i1, 1:j1, end-k2+1:end]
    padded_u[1:i1, end-j2+1:end, end-k2+1:end] = u[1:i1, end-j2+1:end, end-k2+1:end]
    padded_u[end-i2+1:end, 1:j1, end-k2+1:end] = u[end-i2+1:end, 1:j1, end-k2+1:end]
    padded_u[end-i2+1:end, end-j2+1:end, end-k2+1:end] = u[end-i2+1:end, end-j2+1:end, end-k2+1:end]
end

function center_zeropad3(
    u::AbstractArray{T,3}, padded_shape::NTuple{3,Int}) where T

    padded_u = zeros(eltype(u), padded_shape)
    center_zeropad3!(padded_u, u)

    return padded_u
end

norm(x::Vector{AbstractArray{T,N}}) where T where N = sqrt(sum(sum(abs2, ele) for ele in x))

function copy!(dest::Vector{T}, src::Vector{T}) where T<:AbstractArray{S,N} where S where N
    @assert length(dest) == length(src)
    for (d,s) in zip(dest, src)
        copy!(d, s)
    end
end

function block_soft_threshold(x::T, norm_x::T, κ::T)::T where T<:AbstractFloat
    return x * max(T(1.0) - T(κ) / norm_x, T(0))
end


text:
using Documenter, LombScargle

makedocs(
    modules = [LombScargle],
    format = :html,
    sitename = "LombScargle",
)

deploydocs(
    repo = "github.com/JuliaAstro/LombScargle.jl.git",
    target = "build",
    deps = nothing,
    make = nothing,
)


text:
using QuantumLab
using Base.Test

# write your own tests here
@test 1 == 1

# test readGeometryXYZ
@test readGeometryXYZ("h2o.xyz").atoms[2].element.symbol == "O"
@test -1.42 < readGeometryXYZ("h2o.xyz").atoms[3].position.x < -1.41

# test AtomModule
@test Element("C") == Element("C")

# test BaseModule
mqn = MQuantumNumber(0,0,0)
for (mqn in MQuantumNumbers(LQuantumNumber("D"))) end
@test mqn == MQuantumNumber(2,0,0)
@test distance(Position(0,0,1),Position(1,0,0)) == sqrt(2)

# test readBasisSetTX93
@test readBasisSetTX93("STO-3G.tx93").definitions[Element("C")][1].lQuantumNumber.symbol == "S"
@test readBasisSetTX93("STO-3G.tx93").definitions[Element("C")][1].primitives[1].exponent == 71.616837

# test computeBasis
bas = computeBasis(readBasisSetTX93("STO-3G.tx93"),readGeometryXYZ("h2o.xyz"))
@test -0.74 < bas.contractedBFs[3].primitiveBFs[1].center.y < -0.73

# test IntegralsModule
normalize!(bas)
@test_approx_eq 0.0 computeMatrixOverlap(bas)[4,1]
@test_approx_eq_eps 29.00 computeMatrixKinetic(bas)[2,2] 1e-2
@test_approx_eq_eps 0.16175843 IntegralsModule.FIntegral(2,0.3) 1e-8
@test_approx_eq -π IntegralsModule.NuclearAttractionIntegral(PrimitiveGaussianBasisFunction(Position(0,0,0),1.,MQuantumNumber(1,0,0)),PrimitiveGaussianBasisFunction(Position(0,0,0),1.,MQuantumNumber(1,0,0)),Atom(Element("C"),Position(0,0,0)))

max_stars_repo_name: tkelman/QuantumLab

max_issues_repo_name: tkelman/QuantumLab

max_forks_repo_name: tkelman/QuantumLab


text:
using Test
using ODEFilters

using TimerOutputs
const to = TimerOutput()
macro timedtestset(name, code)
    return esc(:(@timeit to $name @testset $name $code))
end

@testset "ODEFilters" begin

    println("Correctness")
    @timedtestset "Correctness" begin include("correctness.jl") end

    println("Filtering")
    @timedtestset "Filtering" begin include("filtering.jl") end

    println("Convergence")
    @timedtestset "Convergence" begin include("convergence.jl") end

    println("Priors")
    @timedtestset "Priors" begin include("priors.jl") end

    println("Solution")
    @timedtestset "Solution" begin include("solution.jl") end

    println("Diffusions")
    @timedtestset "Diffusions" begin include("diffusions.jl") end

    println("State Initialization")
    @timedtestset "State Initialization" begin include("state_init.jl") end

    println("Preconditioning")
    @timedtestset "Preconditioning" begin include("preconditioning.jl") end

    println("Smoothing")
    @timedtestset "Smoothing" begin include("smoothing.jl") end

    println("Errors")
    @timedtestset "Errors" begin include("errors.jl") end

    println("IEKS")
    @timedtestset "IEKS" begin include("ieks.jl") end

    println("Specific Problems")
    @timedtestset "Specific Problems" begin include("specific_problems.jl") end
end

println("\nTimer Summary:")
display(to)
println("\n")


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module StirlingNumbers

using Nemo, IterTools, Triangles

export ModuleStirlingNumbers
export T132393, L132393, V132393, M132393
export T048993, L048993, V048993, M048993
export T271703, L271703, V271703, M271703
export T094587, L094587, V094587, M094587
export T008279, L008279, V008279, M008279
# export Lah, StirlingSet, StirlingCycle

"""

* T132393, L132393, V132393, M132393, T048993, L048993, V048993, M048993, T271703, L271703, V271703, M271703, T094587, L094587, V094587, M094587, T008279, L008279, V008279, M008279
"""
const ModuleStirlingNumbers = ""

"""

Recurrence for A132393, StirlingCycle numbers.
"""
function R132393(n::Int, k::Int, prevrow::Function)
    (k == 0 && n == 0) && return fmpz(1)
    (n - 1) * prevrow(k) + prevrow(k - 1)
end
"""

Recurrence for A048993, StirlingSet numbers.
"""
function R048993(n::Int, k::Int, prevrow::Function)
    (k == 0 && n == 0) && return fmpz(1)
    k * prevrow(k) + prevrow(k - 1)
end
"""

Recurrence for A271703, Lah numbers.
"""
function R271703(n::Int, k::Int, prevrow::Function)
    (k == 0 && n == 0) && return fmpz(1)
    (k - 1 + n) * prevrow(k) + prevrow(k - 1)
end
"""

Recurrence for A094587, (Rising factorials).
"""
function R094587(n::Int, k::Int, prevrow::Function)
    (k == 0 && n == 0) && return fmpz(1)
    (n - k) * prevrow(k) + prevrow(k - 1)
end
"""

Recurrence for A008279. Number of permutations of n things k at a time. (Falling factorials)
"""
function R008279(n::Int, k::Int, prevrow::Function)
    (k == 0 && n == 0) && return fmpz(1)
    prevrow(k) + k * prevrow(k - 1)
end
"""

Iterates over the first n rows of `A132393`.
"""
T132393(n) = RecTriangle(n, R132393)
"""

Iterates over the first n rows of `A048993`.
"""
T048993(n) = RecTriangle(n, R048993)
"""

Iterates over the first n rows of `A271703`.
"""
T271703(n) = RecTriangle(n, R271703)
"""

Iterates over the first n rows of `A094587`.
"""
T094587(n) = RecTriangle(n, R094587)
"""

Iterates over the first n rows of `A008279`.
"""
T008279(n) = RecTriangle(n, R008279)
"""

Lists the first n rows of `A132393` by concatenating.
"""
L132393(n) = vcat(T132393(n)...)
"""

Lists the first n rows of `A048993` by concatenating.
"""
L048993(n) = vcat(T048993(n)...)
"""

Lists the first n rows of `A271703` by concatenating.
"""
L271703(n) = vcat(T271703(n)...)
"""

Lists the first n rows of `A094587` by concatenating.
"""
L094587(n) = vcat(T094587(n)...)
"""

Lists the first n rows of `A008279` by concatenating.
"""
L008279(n) = vcat(T008279(n)...)
"""

Return the triangular array as a square matrix with dim rows.
"""
M132393(dim) = fromΔ(T132393(dim))
"""

Return the triangular array as a square matrix with dim rows.
"""
M048993(dim) = fromΔ(T048993(dim))
"""

Return the triangular array as a square matrix with dim rows.
"""
M271703(dim) = fromΔ(T271703(dim))
"""

Return the triangular array as a square matrix with dim rows.
"""
M094587(dim) = fromΔ(T094587(dim))
"""

Return the triangular array as a square matrix with dim rows.
"""
M008279(dim) = fromΔ(T008279(dim))
"""

Return row n of A132393 based on the iteration `T132393`(n).
"""
V132393(n) = nth(T132393(n + 1), n + 1)
"""

Return row n of A048993 based on the iteration `T048993`(n).
"""
V048993(n) = nth(T048993(n + 1), n + 1)
"""

Return row n of A271703 based on the iteration `T271703`(n).
"""
V271703(n) = nth(T271703(n + 1), n + 1)
"""

Return row n of A094587 based on the iteration `T094587`(n).
"""
V094587(n) = nth(T094587(n + 1), n + 1)
"""

Return row n of A008279 based on the iteration `T008279`(n).
"""
V008279(n) = nth(T008279(n + 1), n + 1)

#START-TEST-########################################################

using Test, SeqUtils

function test()

    # From JuliaMath/Combinatorics/numbers.jl
    function stirlings1(n::Int, k::Int, signed::Bool=false)
        signed == true && return (-1)^(n - k) * stirlings1(n, k)
        if n < 0 throw(DomainError(n, "n must be nonnegative"))
        elseif n == k == 0 return 1
        elseif n == 0 || k == 0  return 0
        elseif n == k return 1
        elseif k == 1 return factorial(n-1)
        elseif k == n - 1 return binomial(n, 2)
        elseif k == n - 2 return div((3 * n - 1) * binomial(n, 3), 4)
        elseif k == n - 3 return binomial(n, 2) * binomial(n, 4)
        end
        return (n - 1) * stirlings1(n - 1, k) + stirlings1(n - 1, k - 1)
    end

    # From JuliaMath/Combinatorics/numbers.jl
    function stirlings2(n::Int, k::Int)
        if n < 0 throw(DomainError(n, "n must be nonnegative"))
        elseif n == k == 0 return 1
        elseif n == 0 || k == 0 return 0
        elseif k == n - 1 return binomial(n, 2)
        elseif k == 2  return 2^(n-1) - 1
        end
        return k * stirlings2(n - 1, k) + stirlings2(n - 1, k - 1)
    end

    @testset "Stirling1" begin

        for (n, S1) ∈ enumerate(T132393(10))
            SC = [stirlings1(n - 1, k) for k ∈ 0:n-1]
            @test S1 == SC
        end
    end

    @testset "Stirling2" begin

        for (n, S2) ∈ enumerate(T048993(10))
            SC = [stirlings2(n - 1, k) for k ∈ 0:n-1]
            @test S2 == SC
        end
    end
end # test

function demo()
    for row ∈ T132393(9)
        Println(row)
    end
    println()

    for n ∈ 0:8
        print(n, ": "); Println(V132393(n))
    end
    println()

    reduce(vcat, [row for row in T132393(6)]) |> Println
    [row for row in T132393(6)] |> println
    [row for row in T132393(6)] |> Println
    for row ∈ T132393(6) Println(row) end

    ShowAsΔ(T132393(9), "\t")

    Println(L132393(9))
    println()

    ShowAsMatrix(T132393(9))
    println()

    println(M132393(9))
    println()

    println(V132393(8))
    println()

    # ===================================

    for row ∈ T048993(9)
        Println(row)
    end
    println()

    for n ∈ 0:8
        print(n, ": "); Println(V048993(n))
    end
    println()

    ShowAsΔ(T048993(9), "\t")

    Println(L048993(9))
    println()

    ShowAsMatrix(T048993(9))
    println()

    println(M048993(9))
    println()

    println(V048993(8))
    println()

    # ===================================

    for row ∈ T008279(9)
        Println(row)
    end
    println()

    for n ∈ 0:8
        print(n, ": "); Println(V008279(n))
    end
    println()

    ShowAsΔ(T008279(9), "\t")

    println(L008279(9))
    println()

    ShowAsMatrix(T008279(9))
    println()

    println(M008279(9))
    println()

    println(V008279(8))
    println()
end

"""

(for row ∈ T132393(500) end)
    0.246288 seconds (379.25 k allocations: 8.748 MiB)
"""
function perf()
    GC.gc()
    @time (for row ∈ T132393(500)
    end)
    GC.gc()
    @time (for row ∈ T048993(500)
    end)
    GC.gc()
    @time (for row ∈ T271703(500)
    end)
end

function main()
    test()
    demo()
    perf()
end

main()

end # module

#=
1;
0,    1;
0,    1,     1;
0,    2,     3,     1;
0,    6,    11,     6,    1;
0,   24,    50,    35,   10,    1;
0,  120,   274,   225,   85,   15,   1;
0,  720,  1764,  1624,  735,  175,  21,  1;
0, 5040, 13068, 13132, 6769, 1960, 322, 28, 1;

1
0  1
0  1   1
0  1   3     1
0  1   7     6     1
0  1   15    25    10     1
0  1   31    90    65     15     1
0  1   63    301   350    140    21    1
0  1   127   966   1701   1050   266   28   1

1
0  1
0  2       1
0  6       6       1
0  24      36      12      1
0  120     240     120     20      1
0  720     1800    1200    300     30      1
0  5040    15120   12600   4200    630     42    1
0  40320   141120  141120  58800   11760   1176  56   1

1
1       1
2       2       1
6       6       3       1
24      24      12      4       1
120     120     60      20      5      1
720     720     360     120     30     6     1
5040    5040    2520    840     210    42    7    1
40320   40320   20160   6720    1680   336   56   8    1

1
1  1
1  2  2
1  3  6   6
1  4  12  24   24
1  5  20  60   120  120
1  6  30  120  360  720   720
1  7  42  210  840  2520  5040  5040
1  8  56  336  1680 6720  20160 40320  40320
=#

max_stars_repo_name: OpenLibMathSeq/Sequences

max_issues_repo_name: OpenLibMathSeq/Sequences


text:
using PyPlot

"""
    plotSpectrum(X, fin, col, nbin=8, n=3)

Plot a smoothed spectrum
"""
function plotSpectrum(X, fin, col, nbin=8, n=3)
    f, p = logsmooth(X, fin, nbin, n)
    semilogx(f, p, color=col, linewidth=0.5)
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
const _CURR_WORKER = GitWorker[]

function gw_curr()
    isempty(_CURR_WORKER) && error("No current GitWorker set (see `gw_curr(gw::GitWorker)`)")
    return first(_CURR_WORKER)
end

function gw_curr(gw::GitWorker)
    empty!(_CURR_WORKER)
    push!(_CURR_WORKER, gw)
    return gw
end

max_stars_repo_name: josePereiro/GitWorkers

max_issues_repo_name: josePereiro/GitWorkers

max_forks_repo_name: josePereiro/GitWorkers


text:
let
    eta = 0.9

    function f1(x)
        (1.0 / 2.0) * (x[1]^2 + eta * x[2]^2)
    end

    function g1(x, storage)
        storage[1] = x[1]
        storage[2] = eta * x[2]
    end

    function h1(x, storage)
        storage[1, 1] = 1.0
        storage[1, 2] = 0.0
        storage[2, 1] = 0.0
        storage[2, 2] = eta
    end

    results = optimize(f1, g1, h1, [127.0, 921.0])
    @assert Optim.g_converged(results)
    @assert norm(Optim.minimizer(results) - [0.0, 0.0]) < 0.01

    results = optimize(f1, g1, [127.0, 921.0])
    @assert Optim.g_converged(results)
    @assert norm(Optim.minimizer(results) - [0.0, 0.0]) < 0.01

    results = optimize(f1, [127.0, 921.0])
    @assert Optim.g_converged(results)
    @assert norm(Optim.minimizer(results) - [0.0, 0.0]) < 0.01

    results = optimize(f1, [127.0, 921.0], autodiff = true)
    @assert Optim.g_converged(results)
    @assert norm(Optim.minimizer(results) - [0.0, 0.0]) < 0.01

    # tests for bfgs_initial_invH
    initial_invH = zeros(2,2)
    h1([127.0, 921.0],initial_invH)
    initial_invH = diagm(diag(initial_invH))
    results = optimize(DifferentiableFunction(f1, g1), [127.0, 921.0], BFGS(), OptimizationOptions(),
                       initial_invH = initial_invH)
    @assert Optim.g_converged(results)
    @assert norm(Optim.minimizer(results) - [0.0, 0.0]) < 0.01
end


text:
using CSV

cd("/Users/kiristern/Documents/GitHub/SDM/data/pred_prey")
df = CSV.read("Mustela_erminea.csv")

#View each sp
# show(unique(df, :species), allrows=true)

#Count number of unique sp observations
show(by(df, :species, :species => length), allrows=true)

max_stars_repo_name: kiristern/SDM

max_issues_repo_name: kiristern/SDM

max_forks_repo_name: kiristern/SDM


text:
# Method of lines discretization scheme

function interface_errors(depvars, indvars, discretization)
    for x in indvars
        @assert findfirst(dxs -> isequal(x, dxs[1].val), discretization.dxs) !== nothing "Variable $x has no step size"
    end
end

function SciMLBase.symbolic_discretize(pdesys::PDESystem, discretization::MethodOfLines.MOLFiniteDifference{G}) where G
    pdeeqs = [eq.lhs - eq.rhs ~ 0 for eq in pdesys.eqs]
    bcs = pdesys.bcs
    domain = pdesys.domain

    t = discretization.time

    depvar_ops = map(x->operation(x.val),pdesys.depvars)
    # Get all dependent variables in the correct type
    alldepvars = get_all_depvars(pdesys, depvar_ops)
    alldepvars = filter(u -> !any(map(x-> x isa Number, arguments(u))), alldepvars)
    # Get all independent variables in the correct type, removing time from the list
    allindvars = remove(collect(filter(x->!(x isa Number), reduce(union, filter(xs->(!isequal(xs, [t])), map(arguments, alldepvars))))), t)
    #@show allindvars, typeof.(allindvars)

    interface_errors(alldepvars, allindvars, discretization)
    # @show alldepvars
    # @show allindvars

    # Get tspan
    tspan = nothing
    # Check that inputs make sense
    if t !== nothing
        tdomain = pdesys.domain[findfirst(d->isequal(t.val, d.variables), pdesys.domain)]
        @assert tdomain.domain isa DomainSets.Interval
        tspan = (DomainSets.infimum(tdomain.domain), DomainSets.supremum(tdomain.domain))
    end
    alleqs = []
    bceqs = []
    # * We wamt to do this in 2 passes
    # * First parse the system and BCs, replacing with DiscreteVariables and DiscreteDerivatives
    # * periodic parameters get type info on whether they are periodic or not, and if they join up to any other parameters
    # * Then we can do the actual discretization by recursively indexing in to the DiscreteVariables

    # Create discretized space and variables
    s = DiscreteSpace(domain, alldepvars, allindvars, discretization)
    # Generate finite difference weights
    derivweights = DifferentialDiscretizer(pdesys, s, discretization)
    # Create a map of each variable to their boundary conditions and get the initial condition
    boundarymap, u0 = BoundaryHandler(pdesys.bcs, s, depvar_ops, tspan, derivweights)
    # Generate a map of each variable to whether it is periodic in a given direction
    pmap = PeriodicMap(boundarymap, s)
    # Get the interior and variable to solve for each equation
    interiormap = InteriorMap(pdeeqs, boundarymap, s)

    # Loop over equations: different space, grid, independent variables etc for each equation
    # a slightly more efficient approach would be to group equations that have the same
    # independent variables
    for pde in pdeeqs
        # Read the dependent variables on both sides of the equation
        depvars_lhs = get_depvars(pde.lhs, depvar_ops)
        depvars_rhs = get_depvars(pde.rhs, depvar_ops)
        depvars = collect(depvars_lhs ∪ depvars_rhs)
        depvars =  filter(u -> !any(map(x-> x isa Number, arguments(u))), depvars)

        # Read the independent variables
        # ignore if the only argument is [t]
        indvars = Set(filter(xs->!isequal(xs, [t]), map(arguments, depvars)))
        # get all parameters in the equation
        allx̄ = Set(filter(!isempty, map(u->filter(x-> t === nothing || !isequal(x, t.val), arguments(u)), depvars)))
        # Handle the case where there are no independent variables apart from time
        if isempty(allx̄)
            rules = varmaps(s, depvars, CartesianIndex(), Dict([]))

            push!(alleqs, substitute(pde.lhs, rules) ~ substitute(pde.rhs, rules))
        else
            # make sure there is only one set of independent variables per equation
            @assert length(allx̄) == 1
            pdex̄ = first(allx̄)
            @assert length(indvars) == 1

            eqvar = interiormap.var[pde]

            args = params(eqvar, s)
            indexmap = Dict([args[i]=>i for i in 1:length(args)])

            # Handle boundary values appearing in the equation by creating functions that map each point on the interior to the correct replacement rule
            # Generate replacement rule gen closures for the boundary values like u(t, 1)
            boundaryvalfuncs = generate_boundary_val_funcs(s, depvars, boundarymap, indexmap, derivweights)

            # Generate the boundary conditions for the correct variable
            for boundary in reduce(vcat, collect(values(boundarymap[operation(eqvar)])))
                generate_bc_eqs!(bceqs, s, boundaryvalfuncs, interiormap, boundary)
            end

            # Set invalid corner points to zero
            generate_corner_eqs!(bceqs, s, interiormap, pde)

            # Generate the equations for the interior points
            pdeeqs = discretize_equation(pde, interiormap.I[pde], eqvar, depvars, s, derivweights, indexmap, boundaryvalfuncs, pmap)

            push!(alleqs,pdeeqs)
        end
    end

    u0 = !isempty(u0) ? reduce(vcat, u0) : u0
    bceqs = reduce(vcat, bceqs)
    alleqs = reduce(vcat, alleqs)
    alldepvarsdisc = unique(reduce(vcat, vec.(values(s.discvars))))

    # Finalize
    defaults = pdesys.ps === nothing || pdesys.ps === SciMLBase.NullParameters() ? u0 : vcat(u0,pdesys.ps)
    ps = pdesys.ps === nothing || pdesys.ps === SciMLBase.NullParameters() ? Num[] : first.(pdesys.ps)
    # Combine PDE equations and BC equations
    try
        if t === nothing
            # At the time of writing, NonlinearProblems require that the system of equations be in this form:
            # 0 ~ ...
            # Thus, before creating a NonlinearSystem we normalize the equations s.t. the lhs is zero.
            eqs = map(eq -> 0 ~ eq.rhs - eq.lhs, vcat(alleqs, unique(bceqs)))
            sys = NonlinearSystem(eqs, vec(reduce(vcat, vec(alldepvarsdisc))), ps, defaults=Dict(defaults),name=pdesys.name)
            return sys, nothing
        else
            # * In the end we have reduced the problem to a system of equations in terms of Dt that can be solved by an ODE solver.

            sys = ODESystem(vcat(alleqs, unique(bceqs)), t, vec(reduce(vcat, vec(alldepvarsdisc))), ps, defaults=Dict(defaults), name=pdesys.name)
            return sys, tspan
        end
    catch e
        println("The system of equations is:")
        println(vcat(alleqs, unique(bceqs)))
        println()
        println("Discretization failed, please post an issue on https://github.com/SciML/MethodOfLines.jl with the failing code and system at low point count.")
        println()
        rethrow(e)
    end
end

function discretize_equation(pde, interior, eqvar, depvars, s, derivweights, indexmap, boundaryvalfuncs, pmap::PeriodicMap{hasperiodic}) where {hasperiodic}
    return vec(map(interior) do II
        boundaryrules = mapreduce(f -> f(II), vcat, boundaryvalfuncs)
        rules = vcat(generate_finite_difference_rules(II, s, depvars, pde, derivweights, pmap, indexmap), boundaryrules, valmaps(s, eqvar, depvars, II, indexmap))
        substitute(pde.lhs,rules) ~ substitute(pde.rhs,rules)
    end)
end

function SciMLBase.discretize(pdesys::PDESystem,discretization::MethodOfLines.MOLFiniteDifference)
    sys, tspan = SciMLBase.symbolic_discretize(pdesys, discretization)
    try
        simpsys = structural_simplify(sys)
        if tspan === nothing
            return prob = NonlinearProblem(simpsys, ones(length(simpsys.states)); discretization.kwargs...)
        else
            return prob = ODEProblem(simpsys, Pair[], tspan; discretization.kwargs...)
        end
    catch e
        error_analysis(sys, e)
    end
end

function get_discrete(pdesys, discretization)
    domain = pdesys.domain

    t = discretization.time

    depvar_ops = map(x->operation(x.val),pdesys.depvars)
    # Get all dependent variables in the correct type
    alldepvars = get_all_depvars(pdesys, depvar_ops)
    alldepvars = filter(u -> !any(map(x-> x isa Number, arguments(u))), alldepvars)
    # Get all independent variables in the correct type, removing time from the list
    allindvars = remove(collect(filter(x->!(x isa Number), reduce(union, filter(xs->(!isequal(xs, [t])), map(arguments, alldepvars))))), t)
    #@show allindvars, typeof.(allindvars)

    interface_errors(alldepvars, allindvars, discretization)
    # @show alldepvars
    # @show allindvars

    # Get tspan
    tspan = nothing
    # Check that inputs make sense
    if t !== nothing
        tdomain = pdesys.domain[findfirst(d->isequal(t.val, d.variables), pdesys.domain)]
        @assert tdomain.domain isa DomainSets.Interval
        tspan = (DomainSets.infimum(tdomain.domain), DomainSets.supremum(tdomain.domain))
    end
    alleqs = []
    bceqs = []

    # Create discretized space and variables
    s = DiscreteSpace(domain, alldepvars, allindvars, discretization)

    return Dict(vcat([Num(x) => s.grid[x] for x in s.x̄], [Num(u) => s.discvars[u] for u in s.ū]) )
end

function ModelingToolkit.ODEFunctionExpr(pdesys::PDESystem,discretization::MethodOfLines.MOLFiniteDifference)
    sys, tspan = SciMLBase.symbolic_discretize(pdesys, discretization)
    try
        if tspan === nothing
            @assert true "Codegen for NonlinearSystems is not yet implemented."
        else
            simpsys = structural_simplify(sys)
            return ODEFunctionExpr(simpsys)
        end
    catch e
        println("The system of equations is:")
        println(sys.eqs)
        println()
        println("Discretization failed, please post an issue on https://github.com/SciML/MethodOfLines.jl with the failing code and system at low point count.")
        println()
        rethrow(e)
    end
end

function generate_code(pdesys::PDESystem,discretization::MethodOfLines.MOLFiniteDifference,filename="code.jl")
    code = ODEFunctionExpr(pdesys, discretization)
    rm(filename)
    open(filename, "a") do io
        println(io, code)
    end
end


text:
using BinaryBuilder, Pkg

name = "LibSSH2"
version = v"1.9.0"

# Collection of sources required to build LibSSH2
sources = [
   "https://github.com/libssh2/libssh2/releases/download/libssh2-$(version)/libssh2-$(version).tar.gz" =>
   "d5fb8bd563305fd1074dda90bd053fb2d29fc4bce048d182f96eaa466dfadafd",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libssh2*/

BUILD_FLAGS=(
    -DCMAKE_BUILD_TYPE=Release
    -DCRYPTO_BACKEND=mbedTLS
    -DBUILD_SHARED_LIBS=ON
    -DBUILD_EXAMPLES=OFF
    -DBUILD_TESTING=OFF
    -DENABLE_ZLIB_COMPRESSION=OFF
    "-DCMAKE_INSTALL_PREFIX=${prefix}"
    "-DCMAKE_TOOLCHAIN_FILE="${CMAKE_TARGET_TOOLCHAIN}""
)

mkdir build && cd build

cmake .. "${BUILD_FLAGS[@]}"
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(;experimental=true)

# The products that we will ensure are always built
products = [
    LibraryProduct("libssh2", :libssh2),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(Pkg.Types.PackageSpec(name="MbedTLS_jll", version=v"2.24.0")),
]

# Note: we explicitly lie about this because we don't have the new
# versioning APIs worked out in BB yet.
version = v"1.9.1"
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
module BloomFilters
	export BloomFilter, add!
	include("bloom-filter.jl")
end


text:
filename = files[end]
states, statenames, units1 = grabstates(filename)
title = "Coarsened " * grabtitle(filename)
grid = getgrid(filename)

filename2 = files[1]
states2, statenames2, units2 = grabstates(filename2)
title2 = "Coarse " * grabtitle(filename2)
coarsesize = size(states2[1])
topo = (Periodic, Bounded, Bounded)
domain = (x=(0, 1e6), y=(0, 1e6), z=(-3000, 0))
newgrid = RegularCartesianGrid(topology=topo, size= coarsesize; domain...)

coarsestates = interpolate(states, grid, newgrid)
scene = visualize(states, states2, statenames = statenames, statenames2 = statenames2, aspect = (1,1, 32/192), statistics = true, title = title, title2 = title2, units1 = units1, units2 = units2)

max_stars_repo_name: sandreza/Mesoscale

max_issues_repo_name: sandreza/Mesoscale

max_forks_repo_name: sandreza/Mesoscale


text:
# generate a random pos def matrix with eigenvalues between 0.1 and 2
function generate_pos_def_matrix(n::Int64,rng)
	X = rand(rng, n, n)
	# any real square matrix can be QP decomposed into a orthogonal matrix and an uppertriangular matrix R
	Q, R = qr(X)
	eigs = rand(rng ,n) .* (2 .- 0.1) .+ 0.1
	X = Q * Matrix(Diagonal(eigs)) * Q'
	X = 0.5 * (X + X')
	return X
end

function is_numerically_pos_sem_def(X, atol)
	X = X ./ 2
	X = X + X'

	F = eigfact(X)
	if size(find( x-> x < -atol, F[:values]), 1) == 0
		return true
	else
		return false
	end
end

function is_numerically_symmetric(X,atol)
	n = size(X, 2)
	for i = 1:n-1, j = i+1:n
		if abs(X[i, j] - X[j, i]) >= atol
			return false
		end
	end
	return true
end

function find_nonsymmetric_component(X)
	for i = 2:size(X, 1), j = 1:(i - 1)
		if abs(X[i, j] - X[j, i]) > 0.0
			return i, j, abs(X[i, j] - X[j, i])
		end
	end
end

function find_different_elements(A, B)
	if size(A) != size(B)
		error("Matrices are not the same size")
	end
	m, n = size(A)
	diff_el = Array[]
	for iii = 1:m, jjj = 1:n
		if A[iii, jjj] != B[iii, jjj]
			push!(diff_el, [iii, jjj])
		end
	end
	return diff_el
end

function duplicate_sparsity_pattern(A)
	m, n = size(A)
	B = zeros(m, n)
	for iii = 1:m, jjj = 1:n
		if A[iii, jjj] != 0
			B[iii, jjj] = 1
		end
	end
	return B

end

function recreate_sparse_matrix(A)
	rowInd = A.rowval
	colPtr = A.colptr
	val = A.nzval

	#compute column indices
	colInd = zeros(Int64, length(rowInd))
	cval = 1
	for iii = 2:length(colPtr)
		currentPtr = colPtr[iii]
		prevPtr = colPtr[iii - 1]
		colInd[prevPtr:currentPtr - 1] = cval
		cval += 1
	end

	# sort rowInd and vals
	p = sortperm(rowInd)
	return sparse(rowInd, colInd, val, size(A, 1), size(A, 2))
end

# Geometric mean from https://github.com/JuliaStats/StatsBase.jl
function gmean(a::AbstractArray{T}) where T<:Real
	s = 0.0
	n = length(a)
	for i in 1:n
		tmp = a[i]
		if tmp < 0.0
			throw(DomainError())
		elseif tmp == 0.0
			return 0.0
		else
			s += log(tmp)
		end
	end
	return exp(s / n)
end


text:
function chambolle_pock(affine_sets::AffineSets, conic_sets::ConicSets, opt)::CPResult#,
    # warm::WarmStart)::CPResult

    # Initialize parameters
    p = Params()
    p.theta = opt.initial_theta
    p.adapt_level = opt.initial_adapt_level
    p.window = opt.convergence_window
    p.beta = opt.initial_beta
    p.time0 = time()
    p.norm_b = norm(affine_sets.b, 2)
    p.norm_h = norm(affine_sets.h, 2)
    p.norm_c = norm(affine_sets.c, 2)
    p.rank_update, p.stop_reason, p.update_cont = 0, 0, 0
    p.stop_reason_string = "Not optimized"
    p.target_rank = 2 * ones(length(conic_sets.sdpcone))
    p.current_rank = 2 * ones(length(conic_sets.sdpcone))
    p.min_eig = zeros(length(conic_sets.sdpcone))
    arc_list = [ARPACKAlloc(Float64, sdp.sq_side) for (idx, sdp) in enumerate(conic_sets.sdpcone)]
    ada_count = 0

    # Print header
    if opt.log_verbose
        print_header_1()
        print_parameters(opt, conic_sets)
        print_constraints(affine_sets)
        if length(conic_sets.socone) + length(conic_sets.sdpcone) > 0
            print_prob_data(conic_sets)
        end
        print_header_2()
    end

    @timeit "Init" begin

        # Scale objective function
        @timeit "normscale alloc" begin
            c_orig, var_ordering = preprocess!(affine_sets, conic_sets)
            A_orig, b_orig = copy(affine_sets.A), copy(affine_sets.b)
            G_orig, h_orig = copy(affine_sets.G), copy(affine_sets.h)
            rhs_orig = vcat(b_orig, h_orig)
        end

        # Diagonal preconditioning
        @timeit "equilibrate" begin
            if opt.equilibration
                M = vcat(affine_sets.A, affine_sets.G)
                UB = maximum(M)
                LB = minimum(M)
                if LB/UB <= opt.equilibration_limit
                    opt.equilibration = false
                end
            end
            if opt.equilibration_force
                opt.equilibration = true
            end
            if opt.equilibration
                @timeit "equilibrate inner" E, D = equilibrate!(M, affine_sets, opt)
                @timeit "equilibrate scaling" begin
                    M = E * M * D
                    affine_sets.A = M[1:affine_sets.p, :]
                    affine_sets.G = M[affine_sets.p + 1:end, :]
                    rhs = E * rhs_orig
                    affine_sets.b = rhs[1:affine_sets.p]
                    affine_sets.h = rhs[affine_sets.p + 1:end]
                    affine_sets.c = D * affine_sets.c
                end
            end
        end

        # Scale the off-diagonal entries associated with p.s.d. matrices by √2
        @timeit "normscale" norm_scaling(affine_sets, conic_sets)

        # Initialization
        pair = PrimalDual(affine_sets)
        a = AuxiliaryData(affine_sets, conic_sets)
        map_socs!(pair.x, conic_sets, a)
        residuals = Residuals(p.window)

        # Diagonal scaling
        M = vcat(affine_sets.A, affine_sets.G)
        Mt = M'

        # Stepsize parameters and linesearch parameters
        if !opt.approx_norm
            @timeit "svd" if minimum(size(M)) >= 2
                try
                    spectral_norm = Arpack.svds(M, nsv=1)[1].S[1]
                catch
                    println("WARNING: Failed to compute spectral norm of M, shifting to Frobenius norm")
                    spectral_norm = norm(M)
                end
            else
                F = LinearAlgebra.svd!(Matrix(M))
                spectral_norm = maximum(F.S)
            end
        else
            spectral_norm = norm(M)
        end

        # Build struct for storing matrices
        mat = Matrices(M, Mt, affine_sets.c)

        # Initial primal and dual steps
        p.primal_step = 1. / spectral_norm
        p.primal_step_old = p.primal_step
        p.dual_step = p.primal_step

        line_search_flag = true
    end

    # Fixed-point loop
    @timeit "CP loop" for k in 1:opt.max_iter

        # Update iterator
        p.iter = k

        # Primal step
        @timeit "primal" primal_step!(pair, a, conic_sets, mat, opt, p, arc_list, p.iter)

        # Linesearch (dual step)
        if line_search_flag
            @timeit "linesearch" linesearch!(pair, a, affine_sets, mat, opt, p)
        else
            @timeit "dual step" dual_step!(pair, a, affine_sets, mat, opt, p)
        end

        # Compute residuals and update old iterates
        @timeit "residual" compute_residual!(residuals, pair, a, p, affine_sets)

        # Compute optimality gap and feasibility error
        @timeit "gap" compute_gap!(residuals, pair, a, affine_sets, p)

        # Print progress
        if opt.log_verbose && mod(k, opt.log_freq) == 0
            print_progress(residuals, p)
        end

        # Check convergence
        p.rank_update += 1
        if residuals.dual_gap <= opt.tol_primal && residuals.equa_feasibility <= opt.tol_primal

            if convergedrank(p, conic_sets, opt) && soc_convergence(a, conic_sets, pair, opt, p)
                p.stop_reason = 1 # Optimal
                p.stop_reason_string = "Optimal solution found"
                if opt.log_verbose
                    print_progress(residuals, p)
                end

                break

            elseif p.rank_update > p.window
                p.update_cont += 1
                if p.update_cont > 0
                    for (idx, sdp) in enumerate(conic_sets.sdpcone)
                        if p.current_rank[idx] + opt.rank_slack >= p.target_rank[idx]
                            if p.min_eig[idx] > opt.tol_psd
                                p.target_rank[idx] = min(2 * p.target_rank[idx], sdp.sq_side)
                                # p.target_rank[idx] = min(1 + p.target_rank[idx], sdp.sq_side)
                            end
                        end
                    end
                    p.rank_update, p.update_cont = 0, 0
                end
            end

        # Check divergence
        elseif k > p.window && residuals.comb_residual[k - p.window] < residuals.comb_residual[k] && p.rank_update > p.window
            p.update_cont += 1
            if p.update_cont > 50
                full_rank_flag = true
                for (idx, sdp) in enumerate(conic_sets.sdpcone)
                    if p.target_rank[idx] < sdp.sq_side
                        full_rank_flag = false
                    end
                    if p.current_rank[idx] + opt.rank_slack >= p.target_rank[idx]
                        if p.min_eig[idx] > opt.tol_psd
                            p.target_rank[idx] = min(2 * p.target_rank[idx], sdp.sq_side)
                            # p.target_rank[idx] = min(1 + p.target_rank[idx], sdp.sq_side)
                        end
                    end
                end
                p.rank_update, p.update_cont = 0, 0
                if full_rank_flag == true
                    p.stop_reason = 4 # Infeasible
                    p.stop_reason_string = "Problem declared infeasible due to lack of improvement"
                    break
                end
            end

        # Adaptive stepsizes
        elseif residuals.primal_residual[k] > opt.tol_primal && residuals.dual_residual[k] < opt.tol_dual && k > p.window
            ada_count += 1
            if ada_count > opt.adapt_window
                ada_count = 0
                if line_search_flag
                    p.beta *= (1. - p.adapt_level)
                    p.primal_step /= sqrt(1. - p.adapt_level)
                else
                    p.primal_step /= (1. - p.adapt_level)
                    p.dual_step *= (1. - p.adapt_level)
                end

                p.adapt_level *= opt.adapt_decay
            end

        elseif residuals.primal_residual[k] < opt.tol_primal && residuals.dual_residual[k] > opt.tol_dual && k > p.window
            ada_count += 1
            if ada_count > opt.adapt_window
                ada_count = 0
                if line_search_flag
                    p.beta /= (1. - p.adapt_level)
                    p.primal_step *= sqrt(1. - p.adapt_level)
                else
                    p.primal_step *= (1. - p.adapt_level)
                    p.dual_step /= (1. - p.adapt_level)
                end

                p.adapt_level *= opt.adapt_decay
            end
        end

        # time_limit stop condition
        if time() - p.time0 > opt.time_limit
            if opt.log_verbose
                print_progress(residuals, p)
            end
            p.stop_reason = 2 # Time limit
            p.stop_reason_string = "Time limit hit, limit: $(opt.time_limit) time: $(time() - p.time0)"
            break
        end

        # max_iter stop condition
        if opt.max_iter == p.iter
            if opt.log_verbose
                print_progress(residuals, p)
            end
            p.stop_reason = 3 # Iteration limit
            p.stop_reason_string = "Iteration limit of $(opt.max_iter) was hit"
        end
    end

    # Remove diag scaling
    cont = 1
    @inbounds for sdp in conic_sets.sdpcone, j in 1:sdp.sq_side, i in 1:j#j:sdp.sq_side
        if i != j
            pair.x[cont] /= sqrt(2.)
        end
        cont += 1
    end

    # Remove equilibrating
    if opt.equilibration
        pair.x = D * pair.x
        pair.y = E * pair.y
    end

    # Equality feasibility error
    equa_error = A_orig * pair.x - b_orig
    # Inequality feasibility error
    slack_ineq = G_orig * pair.x - h_orig

    # Compute results
    time_ = time() - p.time0

    # Print result
    if opt.log_verbose
        print_result(p.stop_reason,
                     time_,
                     residuals,
                     length(p.current_rank) > 0 ? maximum(p.current_rank) : 0)
    end

    # Post processing
    pair.x = pair.x[var_ordering]

    dual_eq = pair.y[1:length(b_orig)]
    dual_in = pair.y[length(b_orig)+1:end]

    return CPResult(p.stop_reason,
                    p.stop_reason_string,
                    pair.x,
                    # pair.y,
                    dual_eq,
                    dual_in,
                    equa_error,
                    slack_ineq,
                    residuals.equa_feasibility,
                    residuals.ineq_feasibility,
                    residuals.prim_obj,
                    residuals.dual_obj,
                    residuals.dual_gap,
                    time_,
                    sum(p.current_rank)
    )
end

function linesearch!(pair::PrimalDual, a::AuxiliaryData, affine_sets::AffineSets, mat::Matrices, opt::Options, p::Params)::Nothing
    cont = 0
    p.primal_step = p.primal_step * sqrt(1. + p.theta)

    for i in 1:opt.max_linsearch_steps
        cont += 1
        p.theta = p.primal_step / p.primal_step_old

        @timeit "linesearch 1" begin
            a.y_half .= pair.y .+ (p.beta * p.primal_step) .* ((1. + p.theta) .* a.Mx .- p.theta .* a.Mx_old)
        end

        @timeit "linesearch 2" begin
            copyto!(a.y_temp, a.y_half)
            box_projection!(a.y_half, affine_sets, p.beta * p.primal_step)
            a.y_temp .-= (p.beta * p.primal_step) .* a.y_half
        end

        @timeit "linesearch 3" mul!(a.Mty, mat.Mt, a.y_temp)

        # In-place norm
        @timeit "linesearch 4" begin
            a.Mty .-= a.Mty_old
            a.y_temp .-= pair.y_old
            y_norm = norm(a.y_temp)
            Mty_norm = norm(a.Mty)
        end

        if sqrt(p.beta) * p.primal_step * Mty_norm <= opt.delta * y_norm
            break
        else
            p.primal_step *= opt.linsearch_decay
        end
    end

    # Reverte in-place norm
    a.Mty .+= a.Mty_old
    a.y_temp .+= pair.y_old

    copyto!(pair.y, a.y_temp)
    p.primal_step_old = p.primal_step
    p.dual_step = p.beta * p.primal_step

    return nothing
end

function dual_step!(pair::PrimalDual, a::AuxiliaryData, affine_sets::AffineSets, mat::Matrices, opt::Options, p::Params)::Nothing

    @timeit "dual step 1" begin
        a.y_half .= pair.y .+ p.dual_step * (2. * a.Mx .- a.Mx_old)
    end

    @timeit "dual step 2" begin
        copyto!(a.y_temp, a.y_half)
        box_projection!(a.y_half, affine_sets, p.dual_step)
        a.y_temp .-= p.dual_step * a.y_half
    end

    @timeit "linesearch 3" mul!(a.Mty, mat.Mt, a.y_temp)

    copyto!(pair.y, a.y_temp)
    p.primal_step_old = p.primal_step

    return nothing
end

function primal_step!(pair::PrimalDual, a::AuxiliaryData, cones::ConicSets, mat::Matrices, opt::Options, p::Params, arc_list, iter::Int64)::Nothing

    pair.x .-= p.primal_step .* (a.Mty .+ mat.c)

    # Projection onto the p.s.d. cone
    if length(cones.sdpcone) >= 1
        @timeit "sdp proj" psd_projection!(pair.x, a, cones, opt, p, arc_list, iter)
    end

    # Projection onto the second order cone
    if length(cones.socone) >= 1
        @timeit "soc proj" soc_projection!(pair.x, a, cones, opt, p)
    end

    @timeit "linesearch -1" mul!(a.Mx, mat.M, pair.x)

    return nothing
end


text:
import 'lib/src/default/math.jl'
import 'lib/src/default/print.jl'
import 'lib/src/default/sort.jl'
import 'lib/src/default/array.jl'

def helper1
def helper2

def main {

    println(flatten([[1, 2], [3, 4], [5, 6]]))

    arr = [1, 5, 3]
    {a, b, arr[1]} = (1, 2, 3)

    print_str('a=')
    println(a)
    print_str('b=')
    println(b)
    print_str('arr[1]=')
    println(arr[1])

    i = 1
    println(i)      ; 1
    println(i++)    ; 1
    println(++i)    ; 3
    println(i)      ; 3
    unsorted = 'floccinaucinihilipilification'
    qsort(unsorted)
    print_strln(unsorted)

    unsorted = [9,8,7,6,5,4,3,2,1]
    unsorted:2 <- 10
    print_str('Unsorted: ')
    println(unsorted)
    qsort(unsorted)
    print_str('Sorted:   ')
    println(unsorted)
}

def main_5 {
    arr = ['abc', 'def', 'ghi', 'jkl', 'mno', 'qrs', 'tuv', 'wxyz']
    print_strln(arr[2])
    println(arr)
}

;def main {
;    if False println(1)
;    else println(0)
;
;    a = [[1, 2], [3, 4], [5, 6]]
;    b = 1
;
;    println(a)
;    println(a[1])
;    println(a[1][1])
;    b = [0, 9, 2, 4]
;    println(b)
;    x <- b:1
;    println(x)
;    3 -> b:2
;    ; not allowed: b:2 <- 3
;    b[0] = 1
;    println(b)
;
;    print(flatten(a), a, 5, 8)
;}

; main
;def main_3 {
;    a = [1, 2, 3, 4, 5]
;    b = [-3, -2, -1, 0]
;
;    a :<= 6
;    a :<= 7
;    a :<= 8
;    0 =>: a
;
;    print(a)
;    ;for i=0, i < |a|, i=i+1
;    ;    print(a[i])
;
;    x <=: a
;    print(x)
;    a :=> x
;    print(x)
;
;    ;for i=0, i < |a|, i=i+1
;    ;    print(a[i])
;
;    print(a, b)
;
;    ;b :>> a
;    ;print(|a|)
;    ;print(|b|)
;
;    ;a[2] = 9
;    ;print(a[2] + 1)
;    ;print([1, 2, 3, 4][2])
;}

def main_2 {
    print(sqrt(2.0))
    print(pow(4,1) as Int )
    print(pow(4,2) as Float )
    print(pow(4,3) as Int )
    print(pow(4,4) as Int )
    print(pow(4,5.551))
}

def main_1 {

    for x = 1, sqr(x) < 50, x++
        print(sqr(x))

    x = 1 ; set x to 1

    while sqr(x) < 50
        print(sqr(x++))

    helper1
    helper2

    if x > 5 exit(1)

    print(5)
}

def helper1 {
    print(sqr(4))
    print(pow(4,1))
    print(pow(4,2))
    print(pow(4,3))
    print(pow(4,4))

    foo = 1.1
    bar = 5f
    print(foo * bar)

    print(sqr(2.5))
}

def helper2 {
    foo = 1
    bar = 2

    foo = foo + bar
    bar = foo

    {baz1, baz2, baz3} = func(bar-3)    ; call func()
    print(baz1)
    print(baz2)
    print(baz3)

    {baz1, baz2, baz3} = (1, 2, 3)
    print(baz1)
    print(baz2)
    print(baz3)

    print(foo + bar)

    fun_if(4)
}

def func(x) {
    print(x)
    print(5)
    100, 200, 300
}

def fun_if(y) {
    if y > 5
        print(-5)
    else if y > 3
        print(-3)
    else
        print(-2)
}

max_stars_repo_name: jeffreymanzione/jeff-vm-lang

max_issues_repo_name: jeffreymanzione/jeff-vm-lang

max_forks_repo_name: jeffreymanzione/jeff-vm-lang


text: @warn "TODO: Implement DIXON3DQ"


text:
#Similar changes to the parallelpiped as were made to the slab on 11.26.2020. This results in new version being >2.9x faster

function DT_refl_paralpip(t, β::Array{Float64,1}, ndet::Float64, nmed::Float64, rd::Array{Float64,1}, rs::Array{Float64,1}, L::Array{Float64,1})

	# rd is location of detector [x,y]
	# rs is location of source [x,y]
	# L is dimensions of parallelpiped [Lx, Ly, Lz]
	n::Float64 = nmed/ndet
    μa::Float64 = β[1]
    μsp::Float64 = β[2]
    D::Float64 = 1/3μsp
	ν::Float64 = 29.9792345/nmed
	xs::UnitRange{Int64} = -10:10

	x::Float64 = rd[1]
	y::Float64 = rd[2]
	xu::Float64 = rs[1]
	yu::Float64 = rs[2]
	lx::Float64 = L[1]
	ly::Float64 = L[2]
	lz::Float64 = L[3]

	x1l = Array{Float64}(undef, length(xs))
	x2l = Array{Float64}(undef, length(xs))
	y1m = Array{Float64}(undef, length(xs))
	y2m = Array{Float64}(undef, length(xs))
	z1n = Array{Float64}(undef, length(xs))
	z2n = Array{Float64}(undef, length(xs))

	Rt1 = Array{Float64}(undef, length(t))
	Rt2 = Array{Float64}(undef, length(t), length(xs))
	Rt3 = Array{Float64}(undef, length(t), length(xs))
	Rt4 = Array{Float64}(undef, length(t), length(xs))

	if n == 1.0
		A= 1.0
	elseif n > 1.0
		A = 504.332889 - 2641.00214n + 5923.699064n^2 - 7376.355814n^3 +
		 5507.53041n^4 - 2463.357945n^5 + 610.956547n^6 - 64.8047n^7
	else
		A = 3.084635 - 6.531194n + 8.357854n^2 - 5.082751n^3
	end

	zo::Float64 = 1/μsp
	zb::Float64 = 2A*D

	x1l = Float64[2l*lx + 4l*zb + xu for l in xs]
	x2l = Float64[2l*lx + (4l-2)*zb - xu for l in xs]
	y1m = Float64[2m*ly + 4m*zb + yu for m in xs]
	y2m = Float64[2m*ly + (4m-2)*zb - yu for m in xs]
	z1n = Float64[2n*lz + 4n*zb + zo for n in xs]
	z2n = Float64[2n*lz + (4n-2)*zb - zo for n in xs]

	Rt1 = @. exp(-μa*ν*t)/(2*(4π*D*ν)^(3/2)*t^(5/2))

	Rt2 = @. exp(-(x-x1l')^2 / (4D*ν*t)) - exp(-(x-x2l')^2 / (4D*ν*t))

	Rt3 = @. exp(-(y-y1m')^2 / (4D*ν*t)) - exp(-(y-y2m')^2 / (4D*ν*t))

	Rt4 = @. z1n'*exp(-z1n'^2 /(4D*ν*t)) - z2n'*exp(-z2n'^2 /(4D*ν*t))

	Rt = Rt1.*sum(Rt2, dims=2).*sum(Rt3, dims=2).*sum(Rt4, dims=2)
	replace!(Rt, NaN =>0)
	return Rt
end

#=
@benchmark DT_refl_paralpip(0:0.01:10, [0.1,10.0], 1.0, 1.0, [2.5,5.0], [5.0,5.0], [10.0,10.0,10.0])
BenchmarkTools.Trial:
  memory estimate:  1.01 MiB
  allocs estimate:  46
  --------------
  minimum time:     1.903 ms (0.00% GC)
  median time:      2.141 ms (0.00% GC)
  mean time:        2.227 ms (2.61% GC)
  maximum time:     5.770 ms (0.00% GC)
  --------------
  samples:          2245
  evals/sample:     1
  =#

  function DT_refl_paralpip1(t, β::Array{Float64,1}, ndet::Float64, nmed::Float64, rd::Array{Float64,1}, rs::Array{Float64,1}, L::Array{Float64,1})

  	# rd is location of detector [x,y]
  	# rs is location of source [x,y]
  	# L is dimensions of parallelpiped [Lx, Ly, Lz]
  		n::Float64 = nmed/ndet
      	μa::Float64 = β[1]
  		μsp::Float64 = β[2]
      	D::Float64 = 1/3μsp
  		ν::Float64 = 29.9792345/nmed
  		xs::UnitRange{Int64} = -10:10

  		x::Float64 = rd[1]
  		y::Float64 = rd[2]
  		xu::Float64 = rs[1]
  		yu::Float64 = rs[2]
  		lx::Float64 = L[1]
  		ly::Float64 = L[2]
  		lz::Float64 = L[3]

		Rt2 = zeros(Float64, length(t))
		Rt3 = zeros(Float64, length(t))
		Rt4 = zeros(Float64, length(t))

  		if n == 1.0
  			A= 1.0
  		elseif n > 1.0
	  		A = 504.332889 - 2641.00214n + 5923.699064n^2 - 7376.355814n^3 +
  		 	5507.53041n^4 - 2463.357945n^5 + 610.956547n^6 - 64.8047n^7
  		else
  			A = 3.084635 - 6.531194n + 8.357854n^2 - 5.082751n^3
  		end

  		zo::Float64 = 1/μsp
  		zb::Float64 = 2A*D

		Rt1 = @. exp(-μa*ν*t)/(2*(4π*D*ν)^(3/2)*t^(5/2))

		Threads.@threads for n in eachindex(t)
			for m in xs
				x1l = 2m.*lx .+ 4m.*zb .+ xu
		  		x2l = 2m.*lx + (4m-2).*zb .- xu
		  		y1m = 2m.*ly .+ 4m.*zb .+ yu
		  		y2m = 2m.*ly .+ (4m.-2).*zb .- yu
		  		z1n = 2m.*lz .+ 4m.*zb .+ zo
		  		z2n = 2m.*lz .+ (4m.-2).*zb .- zo

  				Rt2[n] += exp(-(x-x1l)^2 / (4D*ν*t[n])) - exp(-(x-x2l)^2 / (4D*ν*t[n]))

  				Rt3[n] += exp(-(y-y1m)^2 / (4D*ν*t[n])) - exp(-(y-y2m)^2 / (4D*ν*t[n]))

  				Rt4[n] += z1n*exp(-z1n^2 /(4D*ν*t[n])) - z2n'*exp(-z2n^2 /(4D*ν*t[n]))

			end
		end
  	return replace!(Rt1.*Rt2.*Rt3.*Rt4, NaN => 0)
  end
#=
  @benchmark DT_refl_paralpip1(0:0.01:10, [0.1,10.0], 1.0, 1.0, [2.5,5.0], [5.0,5.0], [10.0,10.0,10.0])
  BenchmarkTools.Trial:
  memory estimate:  43.34 KiB
  allocs estimate:  30
  --------------
  minimum time:     681.021 μs (0.00% GC)
  median time:      733.641 μs (0.00% GC)
  mean time:        990.380 μs (0.41% GC)
  maximum time:     9.611 ms (49.77% GC)
  --------------
  samples:          5046
  evals/sample:     1
=#

max_stars_repo_name: sammzer/LightPropagation

max_issues_repo_name: sammzer/LightPropagation

max_forks_repo_name: sammzer/LightPropagation


text:
## Packages
using ForwardDiff
using UnPack
using Plots
using LinearAlgebra

## Files
include("Tauchen.jl")
include("NewtonRoot.jl")

################################### Model types #########################
mutable struct AiyagariParameters{T <: Real}
    β::T
    σ::T
    A::T
    N::T
    θ::T
    δ::T
end

mutable struct AiyagariModel{T <: Real, I <: Integer}
    r::T
    w::T
    params::AiyagariParameters{T}
    aGrid::Array{T,1} # Policy grid
    aGridl::Array{T,1} # Grid for a with different ϵ stacked up
    na::I
    ϵGrid::Array{T,1} # Labor shock grid
    nϵ::I
    P_ϵ::Array{T,2} # Transition matrix for labor shock process
end

function AiyagariEGM(
    r::T,
    σ::T = 1.0,
    β::T = 0.96,
    A::T = 1.0,
    N::T = 1.0,
    θ::T = 0.33,
    δ::T = 0.05,
    amin::T = 1e-10,
    amax::T = 18.0,
    na::I = 200,
    nϵ::I = 2) where {T <: Real, I <: Integer}

    ########### Parameters
    params = AiyagariParameters(β,σ,A,N,θ,δ);
    w = A * (1 - θ) * (A * θ / (r + δ)) ^ (θ / (1 - θ));

    ########### Policy grid
    aGrid = collect(LinRange(sqrt(amin), sqrt(amax), na));
    aGrid = aGrid.^2;
    aGridl = repeat(aGrid, nϵ)

    ########### Transition
    ϵGrid = [0.1;1.0];
    P_ϵ = [0.9 0.1; 0.1 0.9];

    ########### Initial guess for consumption
    cGuess = r*aGridl + w*repeat(ϵGrid, inner = na)

    return AiyagariModel(r,w,params,aGrid,aGridl,na,ϵGrid,nϵ,P_ϵ), cGuess
end

function interpC(model::AiyagariModel,
                 cbar::AbstractArray,
                 abar::AbstractArray,
                 a_tmr::T,
                 ϵ::T) where {T <: Real}

    @unpack r,w,na = model

    if a_tmr <= abar[1]
        cbar_interp = (1+r)*a_tmr + w*ϵ;
    else
        np = searchsortedlast(abar,a_tmr);

        (np < na) ? np = np :
        np = na-1

        abar_l, abar_h = abar[np], abar[np+1];
        cbar_l, cbar_h = cbar[np], cbar[np+1];
        cbar_interp = cbar_l+(cbar_h-cbar_l)/(abar_h-abar_l)*(a_tmr-abar_l); # y = y0 + (y1-y0)/(x1-x0)*(x-x0)
    end

    return cbar_interp
end

function updateC(model::AiyagariModel,
                 cbar::AbstractArray,
                 abar::AbstractArray,
                 cbar_interp::AbstractArray)

    @unpack r,w,na,nϵ,aGridl,ϵGrid = model

    cbar = reshape(cbar,na,nϵ);
    abar = reshape(abar,na,nϵ);

    for ϵi = 1:nϵ
        for ai = 1:na
            aϵi = (ϵi-1)*na+ai
            cbar_interp[aϵi] = interpC(model,cbar[:,ϵi],abar[:,ϵi],aGridl[aϵi],ϵGrid[ϵi])
        end
    end

    return cbar_interp
end

function EGM(model::AiyagariModel,
             cpol::AbstractArray,
             abar::AbstractArray)

    @unpack r,w,params,na,nϵ,aGridl,P_ϵ,ϵGrid = model
    @unpack β = params

    U(c) = log(c); # Utility function, for now, assume labor is exogenous
    Uc(c) = ForwardDiff.derivative(U, c); # Marginal utility
    Uc_inv(y) = NewtonRoot(c -> Uc(c) - y, 1); # Inverse of marginal utility

    Uc_cj = Uc.(cpol)
    Uc_cj = reshape(Uc_cj, na, nϵ);

    EUc = (P_ϵ*Uc_cj')';
    EUc = reshape(EUc, na*nϵ);
    cbar = Uc_inv.((1+r)*β*EUc);

    # Compute abar
    for ai = 1:na
        for ϵi = 1:nϵ
            aϵi = (ϵi-1)*na+ai
            abar[aϵi] = (cbar[aϵi]+aGridl[aϵi]-w*ϵGrid[ϵi])/(1+r)
        end
    end

    cbar_interp = zeros(length(cbar))
    cbar_interp = updateC(model,cbar,abar,cbar_interp)

    return cbar_interp
end

function SolveEGM(model::AiyagariModel,
                  cCurrent::AbstractArray,
                  abar::AbstractArray,
                  tol = 1e-10)

    @unpack r,w,na,aGridl,nϵ,ϵGrid = model

    for i = 1:10000
        cNext = EGM(model,cCurrent,abar)
        if (i-1) % 50 == 0
            test = maximum(abs.(cCurrent - cNext));
            #println("iteration: $i"," ")
            if test < tol
                #println("Solved in $i iterations")
                break
            end
        end
        cCurrent = copy(cNext);
    end

    apol = zeros(na*nϵ)
    for ai = 1:na
        for ϵi = 1:nϵ
            aϵi = (ϵi-1)*na+ai
            apol[aϵi] = (1+r)*aGridl[aϵi] + w*ϵGrid[ϵi] - cCurrent[aϵi]
        end
    end

    return cCurrent, apol
end

function MakeTransMatEGM(model::AiyagariModel,
                         apol::AbstractArray,
                         Qmat)

    @unpack na,nϵ,aGridl,P_ϵ = model
    apol = reshape(apol,na,nϵ);
    aGridl = reshape(aGridl,na,nϵ);
    for i = 1:na
        for j = 1:nϵ
            aStar = apol[i,j];
            l = searchsortedlast(aGridl[:,j],aStar);

            (l > 0 && l < na) ? l = l :
                (l == na) ? l = na-1 :
                    l = 1

            p = (aStar-aGridl[l,j])/(aGridl[l+1,j]-aGridl[l,j]);
            p = min(max(p,0.0),1.0);

            sj = (j-1)*na;
            for k = 1:nϵ
                sk = (k-1)*na;
                Qmat[sk+l,sj+i] = (1-p)*P_ϵ[k,j];
                Qmat[sk+l+1,sj+i] = p*P_ϵ[k,j];
            end
        end
    end

    return Qmat
end

function StationaryDistributionEGM(Qmat,λCurrent,tol = 1e-10)
    for i = 1:10000
        λnext = Qmat*λCurrent
        if (i-1) % 50 == 0
            test = maximum(abs.(λCurrent - λnext));
            #println("iteration: $i"," ")
            if test < tol
                #println("Solved in $i iterations")
                break
            end
        end
        λCurrent = copy(λnext);
    end

    return λCurrent
end

function equilibriumEGM(model::AiyagariModel,
                        cpol::AbstractArray,
                        abar::AbstractArray,
                        tol = 1e-5,maxr = 50)

    @unpack params,na,nϵ,aGridl,ϵGrid = model
    @unpack β,δ,θ = params

    Ar = 0.0;

    λ_init = zeros(na*nϵ) .+ 1/(na*nϵ);
    ur,lr = 1/β - 1,-δ

    for rit = 1:maxr
        cstar, apol = SolveEGM(model,cpol,abar);

        Qmat = zeros(na*nϵ, na*nϵ);
        Qmat = MakeTransMatEGM(model,apol,Qmat);
        λ = StationaryDistributionEGM(Qmat,λ_init);

        Ar = dot(aGridl, λ);
        Nbar = dot(λ, repeat(ϵGrid, inner = na));
        K = ((model.r+δ)/(θ*Nbar^(1-θ)))^(1/(θ-1));

        println("This is iteration $rit, and r = $(model.r)")

        if (Ar > K)
            ur = model.r;
            model.r = 1.0/2.0*(lr+ur);
        else
            lr = model.r;
            model.r = 1.0/2.0*(lr+ur);
        end
        println("Bond Supply: ",Ar," ","Bond Demand: ",K)

        if abs(Ar - K) < tol
            println("Market clear!")
            return cstar, apol, λ, Ar, K, model.r
            break
        end
    end

    return println("Markets did not clear")
end

##### Test code
# Check for policy functions
r0 = 0.03;
AA0, cpol0 = AiyagariEGM(r0);
abar0 = zeros(length(cpol0));
cpol_EGM, apol_EGM = SolveEGM(AA0,cpol0,abar0);

plot(AA0.aGrid, AA0.aGrid, label = "", color = :black, linestyle = :dash)
plot!(AA0.aGrid, apol_EGM[1:200,:], label = "ϵ = $(round(AA0.ϵGrid[1], digits = 3))")
plot!(AA0.aGrid, apol_EGM[201:400,:], label = "ϵ = $(round(AA0.ϵGrid[2], digits = 3))")

# Check aggregate demand and aggregate supply of capital
r_vals = LinRange(0.005,0.04,20)
Ar_vals = zeros(length(r_vals))
Nbar_vals = zeros(length(r_vals))
K_vals = zeros(length(r_vals))

for i = 1:length(r_vals)
    r0 = r_vals[i];
    AA0, cpol0 = AiyagariEGM(r0);
    abar0 = zeros(length(cpol0));
    cpol_EGM, apol_EGM = SolveEGM(AA0,cpol0,abar0);

    λ_init_test = zeros(AA0.na*AA0.nϵ) .+ 1/(AA0.na*AA0.nϵ);
    Qmat_test = zeros(AA0.na*AA0.nϵ,AA0.na*AA0.nϵ);
    Qmat_test = MakeTransMatEGM(AA0,apol_EGM,Qmat_test);
    λ_test = StationaryDistributionEGM(Qmat_test,λ_init_test)

    Ar_vals[i] = dot(λ_test, AA0.aGridl)
    Nbar_vals[i] = dot(λ_test, repeat(AA0.ϵGrid, inner = AA0.na));
    K_vals[i] = ((r0+AA0.params.δ)/(AA0.params.θ*Nbar_vals[i]^(1-AA0.params.θ)))^(1/(AA0.params.θ-1));
end

plot(Ar_vals,r_vals, label = "supply",lw = 2, alpha = 0.6,xlim = (2, 14), ylim = (0.0, 0.1))
plot!(K_vals,r_vals,label = "demand",lw = 2, alpha = 0.6,xlim = (2, 14), ylim = (0.0, 0.1))