text:
# analysis.jl

"""
    make_dataframes(resultfolder::String,datafolder::String)

"""
function make_dataframes(resultfolder::String,datafolder::String,assemble=true)
    if assemble
        assemble_all(resultfolder)
    end
    files = get_ratesummaryfiles(resultfolder)
    parts = fields.(files)
    models = get_models(parts)
    labels = get_labels(parts)
    df = Vector{Vector}(undef,0)
    for label in labels
        lfiles = files[label .== get_label.(files)]
        dfl = Vector{Tuple}(undef,0)
        for model in models
            mfiles = lfiles[model .== get_model.(lfiles)]
            dfm = Vector{DataFrame}(undef,0)
            for i in eachindex(mfiles)
                push!(dfm,make_dataframe(joinpath(resultfolder,mfiles[i]),datafolder,isfish(mfiles[i])))
            end
            push!(dfl,("Summary_$(label)_$(model).csv",stack_dataframe(dfm)))
        end
        push!(df,dfl)
    end
    return df
end

statfile_from_ratefile(ratefile) = replace(ratefile, "rates_" => "stats_")

function make_dataframe(ratefile::String,datafolder::String,fish::Bool)
    df = read_dataframe(ratefile)
    df2 = read_dataframe(statfile_from_ratefile(ratefile))
    df = leftjoin(df,df2,on = :Gene,makeunique = true)
    filename = splitpath(ratefile)[end]
    parts = fields(filename)
    G = parse(Int,parts.model)
    insertcols!(df, :Model => fill(G,size(df,1)))
    df  = stack_dataframe(df,G,parts.cond)
    add_moments!(df,datafolder,fish)
end

function augment_dataframe(df,resultfolder,datafolder;fishdata=false)
    dfc = copy(df)
    if dfc[1,:Model] > 1
        dfc = add_burstsize(dfc,resultfolder)
    end
    dfc = add_measures(dfc,resultfolder)
    add_modelmoments!(dfc)
end

function make_measure_df(resultfolder::String)
    files = get_measuresummaryfiles(resultfolder)
    df = Vector{DataFrame}(undef,0)
    for file in files
        parts = fields(file)
        df0 = read_dataframe(joinpath(resultfolder,file))
        insertcols!(df0, :Condition => parts.cond)
        push!(df,df0)
    end
    stack_dataframe(df)
end

function add_measures(df,resultfolder::String)
    dm = make_measure_df(resultfolder)
    leftjoin(df,dm[:,[:Gene,:Nalleles,:Condition,:Deviance,:WAIC,:AIC,:Rhat]], on = [:Gene,:Condition])
end

function add_mean!(df::DataFrame,datafolder,fish::Bool)
    # root = string(split(abspath(datafolder),"data")[1])
    m = Vector{Float64}(undef,length(df.Gene))
    i = 1
    for gene in df.Gene
        m[i] = mean_histogram(get_histogram_rna(string(gene),df[i,:Condition],datafolder,fish))
        i += 1
    end
    insertcols!(df, :Expression => m)
end

function add_moments!(df::DataFrame,datafolder,fish::Bool)
    # root = string(split(abspath(datafolder),"data")[1])
    m = Vector{Float64}(undef,length(df.Gene))
    v = similar(m)
    t = similar(m)
    i = 1
    for gene in df.Gene
        m[i] = mean_histogram(get_histogram_rna(string(gene),df[i,:Condition],datafolder,fish))
        v[i] = var_histogram(get_histogram_rna(string(gene),df[i,:Condition],datafolder,fish))
        t[i] = moment_histogram(get_histogram_rna(string(gene),df[i,:Condition],datafolder,fish),3)
        i += 1
    end
    insertcols!(df, :Expression => m, :Variance => v, :ThirdMoment => t)
end

function add_modelmoments!(df::DataFrame)
    m = Vector{Float64}(undef,length(df.Gene))
    v = similar(m)
    i = 1
    for gene in df.Gene
        m[i] = model2_mean(df.Rate01[i],df.Rate10[i],df.Eject[i],df.Decay[i],df.Nalleles[i])
        v[i] = model2_variance(df.Rate01[i],df.Rate10[i],df.Eject[i],df.Decay[i],df.Nalleles[i])
        i += 1
    end
    insertcols!(df, :Model_Expression => m, :Model_Variance => v)
end

add_burstsize(df,resultfolder::String,cols::Vector{Symbol} = [:Gene,:Condition]) = add_burstsize(df,make_burst_df(resultfolder),cols)

add_burstsize(df,db,cols::Vector{Symbol} = [:Gene,:Condition]) = leftjoin(df,db[:,[:BurstSize,:BurstSD,:Gene,:Condition]],on = cols)

function make_burst_df(resultfolder::String)
    files = get_burstsummaryfiles(resultfolder)
    df = Vector{DataFrame}(undef,0)
    for file in files
        parts = fields(file)
        df0 = read_dataframe(joinpath(resultfolder,file))
        insertcols!(df0, :Condition => parts.cond)
        push!(df,df0)
    end
    stack_dataframe(df)
end

# add_time(csvfile::String,timestamp) = CSV.write(csvfile,add_time!(read_dataframe(csvfile),timestamp))

add_time!(df::DataFrame,timestamp) = insertcols!(df, :Time => timestamp)

"""
    stack_dataframe(df,G,cond)
    stack_dataframe(df2::Vector{DataFrame})

"""
stack_dataframe(df,G,cond) = stack_dataframe(separate_dataframe(df,G,cond))

function stack_dataframe(df2::Vector{DataFrame})
    df = df2[1]
    for i in 2:length(df2)
        df = vcat(df,df2[i])
    end
    return df
end

function separate_dataframe(df,G,cond)
    conds = split(cond,"-")
    nsets = length(conds)
    df2 = Vector{DataFrame}(undef,nsets)
    for i in 1:nsets
        df2[i] = df[:,[1; 2*G*(i-1) + 2 : 2*G*i + 1;2*G*nsets+2: end]]
        rename!(x->split(x,"_")[1],df2[i])
        insertcols!(df2[i], :Condition => fill(string(conds[i]),size(df,1)))
    end
    return df2
end

# function residuals(df)
#     R01 = lm(@formula(log(Rate01) ~  Time*log(Decay)+PC1+PC2+PC3+PC4+PC5),df[(df.Winner .> 1) .& (df.Expression .> 0.),:])
#
# end

function make_Zscore_dataframe(df,conditions::Vector)
    dfc = Array{DataFrame}(undef,length(conditions))
    for i in eachindex(dfc)
        dfc[i] = df[df.Condition .== conditions[i],:]
    end
    df = leftjoin(dfc[1],dfc[2], on = [:Gene,:Time],makeunique = true)
    add_Zscore!(df)
    add_MomentChange!(df)
    df[:,[:Gene,:Time,:ZRate01,:ZRate10,:ZEject,:dExpression,:dVariance,:dThirdMoment,:Winner]]
end

function add_MomentChange!(df)
    insertcols!(df, :dExpression => df.Expression_1 .- df.Expression)
    insertcols!(df, :dVariance => df.Variance_1 .- df.Variance)
    insertcols!(df, :dThirdMoment => df.ThirdMoment_1 .- df.ThirdMoment)
end

function add_Zscore!(df)
    insertcols!(df, :ZRate01 => Difference_Zscore.(df.Rate01_1,df.Rate01,df.Rate01SD_1,df.Rate01SD))
    insertcols!(df, :ZRate10 => Difference_Zscore.(df.Rate10_1,df.Rate10,df.Rate10SD_1,df.Rate10SD))
    insertcols!(df, :ZEject => Difference_Zscore.(df.Eject_1,df.Eject,df.EjectSD,df.EjectSD_1))
end

function add_Zscore_class!(df,threshold=2.)
    insertcols!(df, :On => classify_Zscore.(df.ZRate01,threshold))
    insertcols!(df, :Off => classify_Zscore.(df.ZRate10,threshold))
    insertcols!(df, :Eject => classify_Zscore.(df.ZEject,threshold))
    insertcols!(df, :OnOffEject => df.On .* df.Off .* df.Eject)
end

function classify_Zscore(Z,threshold)
    if Z .> threshold
        return "Up"
    elseif Z .< -threshold
        return "Down"
    else
        return "Null"
    end
end

best_AIC(folder::String) = best_measure(folder,:AIC)

best_WAIC(folder::String) = best_measure(folder,:WAIC)

function best_measure(folder::String,measure::Symbol)
    files = get_measuresummaryfiles(folder)
    parts = fields.(files)
    labels = get_labels(parts)
    conds = get_conds(parts)
    df = Vector{Tuple{String,DataFrame}}(undef,0)
    for label in labels
        for cond in conds
            lfiles = files[(label .== get_label.(files)) .& (cond .== get_cond.(files))]
            dm = Vector{DataFrame}(undef,length(lfiles))
            for i in eachindex(lfiles)
                if isfile(joinpath(folder,lfiles[i]))
                    dm[i] = read_dataframe(joinpath(folder,lfiles[i]))
                    insertcols!(dm[i], :Model => fill(get_model(lfiles[i]),size(dm[i],1)))
                end
            end
            # println(best_measure(dm,measure))
            push!(df,("Winners_$(label)_$(cond)_$(string(measure)).csv",best_measure(dm,measure)))
        end
    end
    return df
end

function best_measure(dfs::Vector,measure::Symbol)
    df = DataFrame(:Gene => [], :Winner => [], measure => [])
    ngenes = Int[]
    for d in dfs
        ngenes = push!(ngenes,length(d[:,:Gene]))
    end
    others = setdiff(eachindex(dfs),argmax(ngenes))
    for d in eachrow(dfs[argmax(ngenes)])
        l = d[measure]
        model = d[:Model]
        for k in others
            dc = dfs[k][dfs[k].Gene .== d.Gene,:]
            if ~isempty(dc)
                if dc[1, measure] < l
                    l = dc[1, measure]
                    model = dc[1, :Model]
                end
            end
        end
        push!(df,Dict(:Gene => d.Gene, :Winner => model, measure => l))
    end
    return df
end

"""
    add_pca(df::DataFrame,files::Vector,npcs::Int,conds::Vector)

    make PCAs from files and add to dataframe df

"""

add_pca(df::DataFrame,files::Vector,npcs::Int,conds::Vector) = add_pca(df,make_pca(files,npcs,conds),makeunique=true)

add_pca(df::DataFrame,dpc::DataFrame) = leftjoin(df,dpc,on = [:Condition,:Gene,:Time],makeunique=true)

"""
    make_pca(files::Vector,npcs::Int)

    Compute PCA

"""

function make_pca(files::Vector,npcs::Int,conds::Vector,times::Vector)
    pca = make_pca(files[1],npcs,conds[1],times[1])
    for i in 2:length(files)
        pca = [pca;make_pca(files[i],npcs,conds[i],times[i])]
    end
    return pca
end

function make_pca(files::Vector,npcs::Int,conds::Vector)
    pca = make_pca(files[1],npcs,conds[1])
    for i in 2:length(files)
        pca = [pca;make_pca(files[i],npcs,conds[i])]
    end
    return pca
end

function make_pca(file::String,npcs::Int,cond)
    r,h = readdlm(file,header=true)
    df = make_pca(r,npcs)
    insertcols!(df, :Condition => cond)
end

function make_pca(file::String,npcs::Int,cond,time)
    r,h = readdlm(file,header=true)
    df = make_pca(r,npcs)
    insertcols!(df, :Condition => cond, :Time => time)
end

function make_pca(m::Matrix,npcs::Int)
    M = fit(PCA,float.(m[:,2:end]),maxoutdim = npcs)
    make_pca(M,m[:,1])
end

function make_pca(M::PCA,genes)
    P = projection(M)
    df = DataFrame(Gene = genes)
    i = 1
    for p in eachcol(P)
        insertcols!(df, Symbol("PC$i") => p)
        i += 1
    end
    return df
end

function make_combinedpca(files::Vector,npcs::Int)
    m,h = readdlm(files[1],header=true)
    for i in 2:length(files)
        r,h = readdlm(files[i],header=true)
        m = [m r[:,2:end]]
    end
    make_pca(m,npcs)
end

function make_dataframe_transient(folder::String,winners::String = "")
    rs = Array{Any,2}(undef,0,8)
    files =readdir(folder)
    n = 0
    for file in files
        if occursin("rate",file)
            t = parse(Float64,split(split(file,"T")[2],"_")[1])
            r,head = readdlm(joinpath(folder,file),',',header=true)
            r = [vcat(r[:,[1,2,3,4,5,10]], r[:,[1,6,7,8,9,10]])  [zeros(size(r)[1]); ones(size(r)[1])]  t*ones(2*size(r)[1])/60.]
            rs = vcat(rs,r)
            n += 1
        end
    end
    if winners != ""
        w = get_winners(winners,2*n)
        return DataFrame(Gene = rs[:,1],on = float.(rs[:,2]),off=float.(rs[:,3]),eject=float.(rs[:,4]),decay=float.(rs[:,5]),yield=float.(rs[:,6]),cond = Int.(rs[:,7]),time = float.(rs[:,8]),winner = w)
    else
        return DataFrame(Gene = rs[:,1],on = float.(rs[:,2]),off=float.(rs[:,3]),eject=float.(rs[:,4]),decay=float.(rs[:,5]),yield=float.(rs[:,6]),cond = Int.(rs[:,7]),time = float.(rs[:,8]))
    end
end

function filter_gene(measurefile,measure,threshold)
    genes = Vector{String}(undef,0)
    measures,header = readdlm(measurefile,',',header=true)
    println(length(measures[:,1]))
    col = findfirst(header[1,:] .== measure)
    for d in eachrow(measures)
        if d[col] > threshold || isnan(d[col])
            push!(genes,d[1])
        end
    end
    println(length(genes))
    return genes
end

function filter_gene_nan(measurefile,measure)
    genes = Vector{String}(undef,0)
    measures,header = readdlm(measurefile,',',header=true)
    println(length(measures[:,1]))
    col = findfirst(header[1,:] .== measure)
    for d in eachrow(measures)
        if isnan(d[col])
            push!(genes,d[1])
        end
    end
    println(length(genes))
    return genes
end
"""
large_rhat(measureile,threshold)

"""
large_rhat(measurefile,threshold) = filter_gene(measurefile,"Rhat",threshold)
"""
    large_deviance(measurefile,threshold)

returns list of genes that have deviance greater than `threshold' in 'measurefile'
"""
large_deviance(measurefile,threshold) = filter_gene(measurefile,"Deviance",threshold)

function deviance(r,cond,n,datafolder,root)
    h,hd = histograms(r,cell,cond,n,datafolder,root)
    deviance(h,hd)
end

function compute_deviance(outfile,ratefile::String,cond,n,datafolder,root)
    f = open(outfile,"w")
    rates,head = readdlm(ratefile,',',header=true)
    for r in eachrow(rates)
        d=deviance(r,cond,n,datafolder,root)
        writedlm(f,[gene d],',')
    end
    close(f)
end

frequency(ron,sd,rdecay) = (ron/rdecay, sd/rdecay)

function burstsize(reject::Float64,roff,covee,covoo,coveo::Float64)
        v = var_ratio(reject,roff,covee,covoo,coveo)
        return reject/roff, sqrt(v)
end

function ratestats(gene,G,folder,cond)
    filestats=joinpath(folder,getfile("param-stats",gene,G,folder,cond)[1])
    filerates = joinpath(folder,getratefile(gene,G,folder,cond)[1])
    rates = readrates(filerates)
    cov = read_covparam(filestats)
    # mu = readmean(filestats)
    # sd = readsd(filestats)
    return rates, cov
end

meanofftime(gene::String,infile,n,method,root) = sum(1 .- offtime(gene,infile,n,method,root))

function meanofftime(r::Vector,n::Int,method::Int)
    if n == 1
        return 1/r[1]
    else
        return sum(1 .- offtime(r,n,method))
    end
end

function offtime(r::Vector,n::Int,method::Int)
    _,_,TI = mat_G_DT(r,n)
    vals,_ = eig_decompose(TI)
    minval = min(minimum(abs.(vals[vals.!=0])),.2)
    offtimeCDF(collect(1.:5/minval),r,n,TI,method)
end

function offtime(gene::String,infile,n,method,root)
    contents,head = readdlm(infile,',',header=true)
    r = float.(contents[gene .== contents[:,1],2:end-1])[1,:]
    offtime(r,n,method)

end

function join_files(file1::String,file2::String,outfile::String,addlabel::Bool=true)
    contents1,head1 = readdlm(file1,',',header=true)   # model G=2
    contents2,head2 = readdlm(file2,',',header=true)   # model G=3
    f = open(outfile,"w")
    if addlabel
        header = vcat(String.(head1[2:end]) .* "_G2",String.(head2[2:end]) .* "_G3")
    else
        header = vcat(String.(head1[2:end]),String.(head2[2:end]))
    end
    header = reshape(permutedims(header),(1,length(head1)+length(head2)-2))
    header = hcat(head1[1],header)
    println(header)
    writedlm(f,header,',')
    for row in 1:size(contents1,1)
        if contents1[row,1] == contents2[row,1]
            contents = hcat(contents1[row:row,2:end],contents2[row:row,2:end])
            contents = reshape(permutedims(contents),(1,size(contents1,2)+size(contents2,2)-2))
            contents = hcat(contents1[row,1],contents)
            writedlm(f,contents,',')
        end
    end
    close(f)
end

function join_files(models::Array,files::Array,outfile::String,addlabel::Bool=true)
    m = length(files)
    contents = Array{Array,1}(undef,m)
    headers = Array{Array,1}(undef,m)
    len = 0
    for i in 1:m
        contents[i],headers[i] = readdlm(files[i],',',header=true)
        len += length(headers[i][2:end])
    end
    f = open(outfile,"w")
    header = Array{String,1}(undef,0)
    for i in 1:m
        if addlabel
            header = vcat(header,String.(headers[i][2:end]) .* "_G$(models[i])")
        else
            header = vcat(header,String.(headers[i][2:end]))
        end
    end
    header = reshape(permutedims(header),(1,len))
    header = hcat(headers[1][1],header)
    println(header)
    writedlm(f,header,',')
    for row in 1:size(contents[1],1)
        content = contents[1][row:row,2:end]
        for i in 1:m-1
            if contents[i][row,1] == contents[i+1][row,1]
                content = hcat(content,contents[i+1][row:row,2:end])
                # content = reshape(permutedims(content),(1,len))
            end
        end
        content = hcat(contents[1][row:row,1],content)
        writedlm(f,[content],',')
    end
    close(f)
end

function sample_non1_genes(infile,n)
    contents,head = readdlm(infile,',',header=true)
    list = Array{String,1}(undef,0)
    for c in eachrow(contents)
        if c[5] != 1
            push!(list,c[1])
        end
    end
    a = StatsBase.sample(list,n,replace=false)
end

function add_best_burst(filein,fileout,filemodel2,filemodel3)
    contents,head = readdlm(filein,',',header=true)
    burst2,head2 = readdlm(filemodel2,',',header=true)
    burst3,head3 = readdlm(filemodel3,',',header=true)
    f = open(fileout,"w")
    head = hcat(head,["mean off period" "bust size"])
    writedlm(f,head,',')
    for row in eachrow(contents)
        if Int(row[end]) == 2
            writedlm(f, hcat(permutedims(row),permutedims(burst2[findfirst(burst2[:,1] .== row[1]),2:3])),',')
        else
            writedlm(f, hcat(permutedims(row),permutedims(burst3[findfirst(burst3[:,1] .== row[1]),2:3])),',')
        end
    end
    close(f)
end

function add_best_occupancy(filein,fileout,filemodel2,filemodel3)
    contents,head = readdlm(filein,',',header=true)
    occupancy2,head2 = readdlm(filemodel2,',',header=true)
    occupancy3,head3 = readdlm(filemodel3,',',header=true)
    f = open(fileout,"w")
    head = hcat(head,["Off -2" "Off -1" "On State" ])
    writedlm(f,head,',')
    for row in eachrow(contents)
        if Int(row[end-2]) == 2
            writedlm(f, hcat(permutedims(row),hcat("NA",permutedims(occupancy2[findfirst(occupancy2[:,1] .== row[1]),2:end]))),',')
        else
            writedlm(f, hcat(permutedims(row),permutedims(occupancy3[findfirst(occupancy3[:,1] .== row[1]),2:end])),',')
        end
    end
    close(f)
end

function prune_file(list,file,outfile,header=true)
    contents,head = readdlm(file,',',header=header)
    f = open(outfile,"w")
    for c in eachrow(contents)
        if c[1] in list
            writedlm(f,[c],',')
        end
    end
    close(f)
end

function replace_yield(G,folder1,folder2,cond1,cond2,outfolder)
    if typeof(G) <: Number
        G = string(G)
    end
    if ~isdir(outfolder)
        mkpath(outfolder)
    end
    files1 = getratefile(folder1,G,cond1)
    files2 = getratefile(folder2,G,cond2)
    for file1 in files1
        gene = get_gene(file1)
        file2 = getratefile(files2,gene)
        outfile = joinpath(outfolder,file2)
        r1 = readrates(joinpath(folder1,file1))
        r2 = readrates(joinpath(folder2,file2))
        r2[end] = r1[end]
        f = open(outfile,"w")
        writedlm(f,[r2],',')
        close(f)
    end

end

"""
assemble_r(G,folder1,folder2,cond1,cond2,outfolder)

Combine rates from two separate fits into a single rate vector

"""

function assemble_r(G,folder1,folder2,cond1,cond2,outfolder)
    if typeof(G) <: Number
        G = string(G)
    end
    if ~isdir(outfolder)
        mkpath(outfolder)
    end
    files1 = getratefile(folder1,G,cond1)
    files2 = getratefile(folder2,G,cond2)
    for file1 in files1
        gene = get_gene(file1)
        file2 = getratefile(files2,gene)
        if file2 != 0
            file2=joinpath(folder2,file2)
        else
            file2=joinpath(folder1,file1)
        end
        name = replace(file1, cond1 => cond1 * "-" * cond2)
        outfile = joinpath(outfolder,name)
        assemble_r(joinpath(folder1,file1),file2,outfile)
    end
end

function  assemble_r(ratefile1,ratefile2,outfile)
    r1 = readrates(ratefile1,2)
    r2 = readrates(ratefile2,2)
    r1[end] = clamp(r1[end],eps(Float64),1-eps(Float64))
    r = vcat(r1[1:end-1],r2[1:end-1],r1[end])
    f = open(outfile,"w")
    writedlm(f,[r],',')
    writedlm(f,[r],',')
    writedlm(f,[r],',')
    writedlm(f,[r],',')
    close(f)
end

function assemble_r(gene,G,folder1,folder2,cond1,cond2,outfolder)
    file1 = getratefile(gene,G,folder1,cond1)[1]
    file2 = getratefile(gene,G,folder2,cond2)[1]
    name = replace(file1, cond1 => cond1 * "-" * cond2)
    println(name)
    outfile = joinpath(outfolder,name)
    println(outfile)
    assemble_r(joinpath(folder1,file1),joinpath(folder2,file2),outfile)
end

function getratefile(files,gene)
    files = files[occursin.("_"*gene*"_",files)]
    if length(files) > 0
        return files[1]
    else
        # println(gene)
        return 0
    end
end

function getratefile(folder,G,cond)
    files = readdir(folder)
    files = files[occursin.("rates_",files)]
    files = files[occursin.("_"*cond*"_",files)]
    files[occursin.("_"*G*"_",files)]
end

getratefile(gene,G,folder,cond) = getfile("rate",gene,G,folder,cond)

function getfile(type,gene::String,G::String,folder,cond)
    files = readdir(folder)
    files = files[occursin.(type,files)]
    files = files[occursin.("_"*gene*"_",files)]
    files = files[occursin.("_"*G*"_",files)]
    files[occursin.("_"*cond*"_",files)]
end

function change_name(folder,oldname,newname)
    files = readdir(folder)
    files = files[occursin.(oldname,files)]
    for file in files
        newfile = replace(file, oldname => newname)
        mv(joinpath(folder,file),joinpath(folder,newfile),force=true)
    end
end

function make_halflife(infile,outfile,col=4)
    f = open(outfile,"w")
    writedlm(f,["Gene" "Halflife"],',')
    contents,rows = readdlm(infile,',',header=true)
    for row = eachrow(contents)
        gene = string(row[1])
        gene = strip(gene,'*')
        h1 = float(row[col])
        h2 = float(row[col+1])
        if h1 > 0 || h2 > 0
            h = (h1 + h2)/(float(h1>0) + float(h2>0))
            writedlm(f,[gene h],',')
        end
    end
    nothing
end

function make_datafiles(infolder,outfolder,label)
    histograms = readdir(infolder)
    if ~isdir(outfolder)
        mkpath(outfolder)
    end
    for file in histograms
        newfile = replace(file,label => "")
        cp(joinpath(infolder,file),joinpath(outfolder,newfile))
    end
end

"""
    histograms(r,cell,cond,n::Int,datafolder,root)

"""
function histograms(rin,cell,cond,G::Int,datafolder,fish,root)
    gene = string(rin[1])
    r = float.(rin[2:end])
    data = data_rna(gene,cond,datafolder,fish,"label",root)
    nalleles = alleles(gene,cell,root)
    model = model_rna(r,[],G,nalleles,.01,[],(),fish)
    likelihoodarray(r,data,model)
end

function get_histogram_rna(gene,cond,datafolder,fish,root)
    if fish
        datapath = FISHpath(gene,cond,datafolder,root)
        h = read_fish(datapath,cond,.98)
    else
        datapath = scRNApath(gene,cond,datafolder,root)
        h = read_scrna(datapath,.99)
    end
    normalize_histogram(h)
end

function get_histogram_rna(gene,cond,datafolder,fish)
    if fish
        datapath = FISHpath(gene,cond,datafolder)
        h = read_fish(datapath,cond,.98)
    else
        datapath = scRNApath(gene,cond,datafolder)
        h = read_scrna(datapath,.99)
    end
    normalize_histogram(h)
end

"""
plot_histogram()
plot_histogram(ratefile::String,datafolder;fish=false,root=".",row=2)

functions to plot data and model predicted histograms

"""

function plot_histogram(ratefile::String,datafolder;fish=false,root=".",row=2)
    r = readrow(ratefile,row)
    println(r)
    parts = fields(ratefile)
    label = parts.label
    cond = parts.cond
    G = parts.model
    data = data_rna(parts.gene,parts.cond,datafolder,fish,parts.label,root)
    model = model_rna(r,r,parse(Int,parts.model),parse(Int,parts.nalleles),.01,[],(),0)
    # model_rna(r,[rateprior[i]],G,nalleles,cv,[fittedparam[i]],fixedeffects,0)
    plot_histogram(data,model)
    return data,model
end

function plot_histogram(gene::String,cell::String,G::Int,cond::String,fish::Bool,ratefile::String,datafolder::String,root::String,yield = -1.)
    rates = readdlm(ratefile,',',header=true)
    r = rates[1][findfirst(rates[1][:,1] .== gene)[1],2:end]
    if yield > 0 && ~fish
        r[end] = yield
    end
    data = data_rna(gene,cond,datafolder,fish,"label",root)
    nalleles = alleles(gene,cell,root)
    model = model_rna(r,[],G,nalleles,.01,[],(),fish,0)
    println(typeof(model))
    println(typeof(data))
    m = plot_histogram(data,model)
    return m,data,model
end

function plot_histogram(gene::String,cell::String,G::String,cond::String,fish::Bool,label::String,ratefolder::String,datafolder::String,nsets::Int,root::String,fittedparam = [1],verbose=false)
    data = data_rna(gene,cond,datafolder,fish,label,root)
    model = model_rna(gene,cell,G,fish,.01,fittedparam,(),label,ratefolder,nsets,root,data,verbose)
    m = plot_histogram(data,model)
    return m,data,model
end

function plot_histogram(data::RNAData{Vector{Int64}, Vector{Array}},model::GMlossmodel)
    h=likelihoodarray(model.rates,data,model)
    for i in eachindex(h)
        figure()
        plot(h[i])
        plot(normalize_histogram(data.histRNA[i]))
        savefig(string(i))
    end
    return h
end
function plot_histogram(data::AbstractRNAData{Array{Array,1}},model)
    h=likelihoodarray(model.rates,data,model)
    figure(data.gene)
    for i in eachindex(h)
        plot(h[i])
        plot(normalize_histogram(data.histRNA[i]))
        savefig(string(i))
    end
    return h
end
function plot_histogram(data::AbstractRNAData{Array{Float64,1}},model)
    h=likelihoodfn(get_param(model),data,model)
    plt = plot(h)
    plot!(plt,normalize_histogram(data.histRNA))
    display(plt)
    return h
end

function plot_histogram(data::RNALiveCellData,model::AbstractGRMmodel)
    h=likelihoodtuple(model.rates,data,model)
    plt1 = plot(h[1])
    plot!(plt1,normalize_histogram(data.OFF))
    plt2 = plot(h[2])
    plot!(plt2,normalize_histogram(data.ON))
    plt3 = plot(h[3])
    plot!(plt3,normalize_histogram(data.histRNA))
    plt = plot(plt1,plt2,plt3,layout = (3,1))
    display(plt)
    return h
end

function plot_histogram(data::TransientRNAData,model::AbstractGMmodel)
    h=likelihoodarray(model.rates,data,model)
    for i in eachindex(h)
        figure(data.gene *":T" * "$(data.time[i])")
        plot(h[i])
        plot(normalize_histogram(data.histRNA[i]))
    end
    return h
end

function plot_histogram(data::RNAData{T1,T2},model::AbstractGMmodel,save = false) where {T1 <: Array, T2 <: Array}
    m=likelihoodarray(model.rates,data,model)
    println("*")
    for i in eachindex(m)
        plt=plot(m[i])
        plot!(normalize_histogram(data.histRNA[i]),show=true)
        if save
            savefig()
        else
            display(plt)
        end
    end
    println(deviance(data,model))
    println(loglikelihood(get_param(model),data,model)[1])
    return m
end

function plot_histogram(data::RNAData,model::AbstractGMmodel)
    h=likelihoodfn(get_param(model),data,model)
    plt=plot(h)
    plot!(normalize_histogram(data.histRNA))
    display(plt)
    return h
end

function plot_model(r,n,nhist,nalleles,yield)
    h= steady_state(r[1:2*n+2],yield,n,nhist,nalleles)
    plt = plot(h)
    display(plt)
    return h
end

function plot_model(r,n,nhist,nalleles)
    h= steady_state(r[1:2*n+2],n,nhist,nalleles)
    plt = plot(h)
    display(plt)
    return h
end


text:
# ---
# title: 1178. Number of Valid Words for Each Puzzle
# id: problem1178
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Hash Table, Bit Manipulation
# link: <https://leetcode.com/problems/number-of-valid-words-for-each-puzzle/description/>
# hidden: true
# ---
#
# With respect to a given `puzzle` string, a `word` is _valid_  if both the
# following conditions are satisfied:
#
#   * `word` contains the first letter of `puzzle`.
#   * For each letter in `word`, that letter is in `puzzle`.
# For example, if the puzzle is "abcdefg", then valid words are "faced",
# "cabbage", and "baggage"; while invalid words are "beefed" (doesn't include
# "a") and "based" (includes "s" which isn't in the puzzle).
#
# Return an array `answer`, where `answer[i]` is the number of words in the
# given word list `words` that are valid with respect to the puzzle
# `puzzles[i]`.
#
#
#
# **Example :**
#
#
#
#     Input:
#     words = ["aaaa","asas","able","ability","actt","actor","access"],
#     puzzles = ["aboveyz","abrodyz","abslute","absoryz","actresz","gaswxyz"]
#     Output: [1,1,3,2,4,0]
#     Explanation:
#     1 valid word for "aboveyz" : "aaaa"
#     1 valid word for "abrodyz" : "aaaa"
#     3 valid words for "abslute" : "aaaa", "asas", "able"
#     2 valid words for "absoryz" : "aaaa", "asas"
#     4 valid words for "actresz" : "aaaa", "asas", "actt", "access"
#     There're no valid words for "gaswxyz" cause none of the words in the list contains letter 'g'.
#
#
#
#
# **Constraints:**
#
#   * `1 <= words.length <= 10^5`
#   * `4 <= words[i].length <= 50`
#   * `1 <= puzzles.length <= 10^4`
#   * `puzzles[i].length == 7`
#   * `words[i][j]`, `puzzles[i][j]` are English lowercase letters.
#   * Each `puzzles[i] `doesn't contain repeated characters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end

ext: jl


text:
type OggDecoder
    sync_state::OggSyncState
    streams::Dict{Clong,OggStreamState}
    pages::Dict{Clong,Vector{OggPage}}
    packets::Dict{Clong,Vector{Vector{UInt8}}}

    function OggDecoder()
        syncref = Ref{OggSyncState}(OggSyncState())
        status = ccall((:ogg_sync_init,libogg), Cint, (Ref{OggSyncState},), syncref)
        dec = new(syncref[], Dict{Clong,OggStreamState}(), Dict{Clong,Vector{OggPage}}(), Dict{Clong,Vector{Vector{UInt8}}}())
        if status != 0
            error("ogg_sync_init() failed: This should never happen")
        end

        # This seems to be causing problems.  :(
        # finalizer(dec, x -> begin
        #     for serial in keys(x.streams)
        #         ogg_stream_destroy(x.streams[serial])
        #     end
        #     ogg_sync_destroy(x.sync_state)
        # end )

        return dec
    end
end

function show(io::IO, x::OggDecoder)
    num_streams = length(x.streams)
    if num_streams != 1
        write(io, "OggDecoder with $num_streams streams")
    else
        write(io, "OggDecoder with 1 stream")
    end
end

function ogg_sync_buffer(dec::OggDecoder, size)
    syncref = Ref{OggSyncState}(dec.sync_state)
    buffer = ccall((:ogg_sync_buffer,libogg), Ptr{UInt8}, (Ref{OggSyncState}, Clong), syncref, size)
    dec.sync_state = syncref[]
    if buffer == C_NULL
        error("ogg_sync_buffer() failed: returned buffer NULL")
    end
    return pointer_to_array(buffer, size)
end

function ogg_sync_wrote(dec::OggDecoder, size)
    syncref = Ref{OggSyncState}(dec.sync_state)
    status = ccall((:ogg_sync_wrote,libogg), Cint, (Ref{OggSyncState}, Clong), syncref, size)
    dec.sync_state = syncref[]
    if status != 0
        error("ogg_sync_wrote() failed: error code $status")
    end
end

function ogg_sync_pageout(dec::OggDecoder)
    syncref = Ref{OggSyncState}(dec.sync_state)
    pageref = Ref{OggPage}(OggPage())
    status = ccall((:ogg_sync_pageout,libogg), Cint, (Ref{OggSyncState}, Ref{OggPage}), syncref, pageref)
    dec.sync_state = syncref[]
    if status == 1
        return pageref[]
    else
        return nothing
    end
end

function ogg_page_serialno(page::OggPage)
    pageref = Ref{OggPage}(page)
    return Clong(ccall((:ogg_page_serialno,libogg), Cint, (Ref{OggPage},), pageref))
end

function ogg_page_eos(page::OggPage)
    pageref = Ref{OggPage}(page)
    return ccall((:ogg_page_eos,libogg), Cint, (Ref{OggPage},), pageref)
end

"""
Send a page in, return the serial number of the stream that we just decoded
"""
function ogg_stream_pagein(dec::OggDecoder, page::OggPage)
    serial = ogg_page_serialno(page)
    if !haskey(dec.streams, serial)
        streamref = Ref{OggStreamState}(OggStreamState())
        status = ccall((:ogg_stream_init,libogg), Cint, (Ref{OggStreamState}, Cint), streamref, serial)
        if status != 0
            error("ogg_stream_init() failed: Unknown failure")
        end
        dec.streams[serial] = streamref[]

        # Also initialize dec.packets and dec.pages for this serial
        dec.pages[serial] = Vector{Vector{UInt8}}()
        dec.packets[serial] = Vector{Vector{UInt8}}()
    end

    # Save the page in dec.pages for posterity
    push!(dec.pages[serial], page)

    streamref = Ref{OggStreamState}(dec.streams[serial])
    pageref = Ref{OggPage}(page)
    status = ccall((:ogg_stream_pagein,libogg), Cint, (Ref{OggStreamState}, Ref{OggPage}), streamref, pageref)
    dec.streams[serial] = streamref[]
    if status != 0
        error("ogg_stream_pagein() failed: Unknown failure")
    end
    return serial
end

function ogg_stream_packetout(dec::OggDecoder, serial::Clong; retry::Bool = false)
    if !haskey(dec.streams, serial)
        return nothing
    end
    streamref = Ref{OggStreamState}(dec.streams[serial])
    packetref = Ref{OggPacket}(OggPacket())
    status = ccall((:ogg_stream_packetout,libogg), Cint, (Ref{OggStreamState}, Ref{OggPacket}), streamref, packetref)
    dec.streams[serial] = streamref[]
    if status == 1
        return packetref[]
    else
        # Is our status -1?  That means we're desynchronized and should try again, at least once
        if status == -1 && !retry
            return ogg_stream_packetout(dec, serial; retry = true)
        end
        return nothing
    end
end

function decode_all_pages(dec::OggDecoder, enc_io::IO; chunk_size::Integer = 4096)
    # Load data in until we have a page to sync out
    while !eof(enc_io)
        page = ogg_sync_pageout(dec)
        while page != nothing
            ogg_stream_pagein(dec, page)
            page = ogg_sync_pageout(dec)
        end

        # Load in up to `chunk_size` of data, unless the stream closes before that
        buffer = ogg_sync_buffer(dec, chunk_size)
        bytes_read = readbytes!(enc_io, buffer, chunk_size)
        ogg_sync_wrote(dec, bytes_read)
    end

    # Do our last pageouts to get the last pages
    page = ogg_sync_pageout(dec)
    while page != nothing
        ogg_stream_pagein(dec, page)
        page = ogg_sync_pageout(dec)
    end
end

"""
File goes in, packets come out
"""
function decode_all_packets(dec::OggDecoder, enc_io::IO)
    # Now, decode all packets for these pages
    for serial in keys(dec.streams)
        packet = ogg_stream_packetout(dec, serial)
        while packet != nothing
            push!(dec.packets[serial], copy(pointer_to_array(packet.packet, packet.bytes)))

            if packet.e_o_s == 1
                delete!(dec.streams, serial)
            end

            packet = ogg_stream_packetout(dec, serial)
        end
    end
end

function load(fio::IO; chunk_size=4096)
    dec = OggDecoder()
    decode_all_pages(dec, fio; chunk_size=chunk_size)
    decode_all_packets(dec, fio)
    return dec.packets
end

function load(file_path::Union{File{format"OGG"},AbstractString}; chunk_size=4096)
    open(file_path) do fio
        return load(fio)
    end
end


text:
using DelimitedFiles

open("C:\\Users\\warner\\Downloads\\Lib80x\\xs_dir.txt") do f
    line = readline(f)
    print(typeof(f))
    print(line)
end

max_stars_repo_name: wamcghee/NautilusMC

max_issues_repo_name: wamcghee/NautilusMC

max_forks_repo_name: wamcghee/NautilusMC


text:
"""
    quadratic_matrix_form(f, vars)

Extract the linear and quadratic matrices of a system of polynomial equations.

### Input

- `f`    -- vector of polynomials
- `vars` -- vector of variables

### Output

Matrices `F1` and `F2` of size `n × n` and `n × n^2` respectively,
such that `F1` contains the linear coefficients of the polynomial vector
field `f` and `F2` the quadratic coefficients, which indexing that respects
the Kronecker power `vars ⊗ vars`.

### Examples

TO-DO
"""
function quadratic_matrix_form(f, vars)
    n = length(vars)
    @assert length(f) == n

    # get the linear part
    F1 = zeros(n, n)
    update!(F1, vars, f)

    # get the quadratic part
    F2 = zeros(n, n^2)
    update!(F2, kron_pow(vars, 2), f)

    return F1, F2
end

function update!(Fk, xk, f)
    n = size(Fk, 1)
    for i in 1:n
        monom = monomials(f[i])
        S = Set()
        for (j, x) in enumerate(xk)
            # skip monomials already seen
            x ∈ S ? continue : push!(S, x)
            idx = findfirst(==(x), monom)
            if !isnothing(idx)
                Fk[i, j] = coefficient(f[i], x)
            end
        end
    end
    return Fk
end


text:
module rocBLAS

using ..ROCArrays

using LinearAlgebra

include("librocblas_types.jl")
include("error.jl")

const _handle = Ref{rocblas_handle}(C_NULL)
function handle()
    if _handle[] == C_NULL
        handle = rocblas_create_handle()
        atexit(()->rocblas_destroy_handle(handle))
        _handle[] = handle
    end
    return _handle[]
end

include("librocblas.jl")
#include("util.jl")
include("wrappers.jl")
include("highlevel.jl")

function version()
    VersionNumber(join(split(rocblas_get_version_string(), '.')[1:3], '.'))
end

end


text:
# tests for the LinearSolver module

function test_linearsolver()

  fname = "input_vals_linear.jl"
  mesh, sbp, eqn, opts = solvePDE(fname)

  @testset "----- Testing Linear Solver -----" begin
    test_dense(mesh, sbp, eqn, opts)
    test_sparsedirect(mesh, sbp, eqn, opts)
    test_petscmat(mesh, sbp, eqn, opts)
    test_petscmatfree(mesh, sbp, eqn, opts)
    test_petscmat_matfree(mesh, sbp, eqn, opts)
  end

end

function test_dense(mesh, sbp, eqn, opts)

  pc = PCNone(mesh, sbp, eqn, opts)
  lo = DenseLO(pc, mesh, sbp, eqn, opts)
  ls = StandardLinearSolver(pc, lo, eqn.comm, opts)

  ctx_residual = (evalResidual,)
  t = 0.0
  calcLinearOperator(ls, mesh, sbp, eqn, opts, ctx_residual, t)
  rand!(lo.A)
  for i=1:mesh.numDof
    lo.A[i, i] = 5  # make non-singular
  end

  A_orig = copy(lo.A)
  b = rand(mesh.numDof)

  # test product
  c = lo.A*b
  c2 = zeros(c)
  applyLinearOperator(ls.lo, mesh, sbp, eqn, opts, ctx_residual, t, b, c2)
  @test isapprox( norm(c - c2), 0.0) atol=1e-12

  # test transpose product
  c = lo.A.'*b
  applyLinearOperatorTranspose(ls.lo, mesh, sbp, eqn, opts, ctx_residual, t, b, c2)
  @test isapprox( norm(c - c2), 0.0) atol=1e-12

  # test solve
  x = lo.A\b
  x2 = zeros(x)
  linearSolve(ls, b, x2)
  @test isapprox( norm(x2 - x), 0.0) atol=1e-12

  # test PC
  x3 = zeros(x2)
  applyPC(ls, mesh, sbp, eqn, opts, t, b, x3)
  @test isapprox( norm(x3 - x2), 0.0) atol=1e-12

  # test transpose solve
  x = A_orig.'\b
  x2 = zeros(x)
  linearSolveTranspose(ls, b, x2)
  @test isapprox( norm(x2 - x), 0.0) atol=1e-12

  # test PC
  x3 = zeros(x2)
  applyPCTranspose(ls, mesh, sbp, eqn, opts, t, b, x3)
  @test isapprox( norm(x3 - x2), 0.0) atol=1e-12

  # test product (again)
  c = A_orig*b
  c2 = zeros(c)
  @test_throws Exception  applyLinearOperator(ls.lo, mesh, sbp, eqn, opts, ctx_residual, t, b, c2)
  #=
  @test isapprox( norm(c - c2), 0.0) atol=1e-12
  println("c = \n", c)
  println("c2 = \n", c2)
  println("diff = \n", c2 - c)
  println("ipiv = \n", ls.lo.ipiv)
  =#

  # transpose product throws
  @test_throws Exception  applyLinearOperatorTranspose(ls.lo, mesh, sbp, eqn, opts, ctx_residual, t, b, c2)

  # test that the factorization is not performed again

  x = A_orig\b
  x3 = zeros(x2)
  linearSolve(ls, b, x3)
  @test isapprox( norm(x3 - x), 0.0) atol=1e-12
  @test ( ls.lo.nfactorizations )== 1
  @test ( ls.lo.nsolves )== 3

  free(ls)

  return nothing
end

function test_sparsedirect(mesh, sbp, eqn, opts)

  # test sparse direct solve
  println("testing sparse direct linear solver")
  pc = PCNone(mesh, sbp, eqn, opts)
  lo = SparseDirectLO(pc, mesh, sbp, eqn, opts)
  ls = StandardLinearSolver(pc, lo, eqn.comm, opts)

  ctx_residual = (evalResidual,)
  t = 0.0
  calcLinearOperator(ls, mesh, sbp, eqn, opts, ctx_residual, t)
  rand!(lo.A.nzval)
  for i=1:mesh.numDof
    lo.A[i, i] = 5  # make non-singular
  end

  A_orig = copy(lo.A)
  b = rand(mesh.numDof)

  # test product
  c = lo.A*b
  c2 = zeros(c)
  applyLinearOperator(ls.lo, mesh, sbp, eqn, opts, ctx_residual, t, b, c2)
  @test isapprox( norm(c - c2), 0.0) atol=1e-12

  # test transpose product
  c = A_orig.'*b
  applyLinearOperatorTranspose(ls.lo, mesh, sbp, eqn, opts, ctx_residual, t, b, c2)
  @test isapprox( norm(c - c2), 0.0) atol=1e-12

  # test solve
  x = factorize(lo.A)\b
  x2 = zeros(x)
  linearSolve(ls, b, x2)
  @test isapprox( norm(x2 - x), 0.0) atol=1e-12

  # test PC
  x3 = zeros(x2)
  applyPC(ls, mesh, sbp, eqn, opts, t, b, x3)
  @test isapprox( norm(x3 - x2), 0.0) atol=1e-12

  # test transpose solve
  x = A_orig.'\b
  x2 = zeros(x)
  linearSolveTranspose(ls, b, x2)
  @test isapprox( norm(x2 - x), 0.0) atol=1e-12

  # test PC
  x3 = zeros(x2)
  applyPCTranspose(ls, mesh, sbp, eqn, opts, t, b, x3)
  @test isapprox( norm(x3 - x2), 0.0) atol=1e-12

  # test product (again)
  c = lo.A*b
  c2 = zeros(c)
  applyLinearOperator(ls.lo, mesh, sbp, eqn, opts, ctx_residual, t, b, c2)
  @test isapprox( norm(c - c2), 0.0) atol=1e-12

  # test transpose product (again)
  c = A_orig.'*b
  applyLinearOperatorTranspose(ls.lo, mesh, sbp, eqn, opts, ctx_residual, t, b, c2)
  @test isapprox( norm(c - c2), 0.0) atol=1e-12

  # test solve (again)
  x = factorize(lo.A)\b
  x2 = zeros(x)
  linearSolve(ls, b, x2)
  @test isapprox( norm(x2 - x), 0.0) atol=1e-12
  @test ( ls.lo.nfactorizations )== 1
  @test ( ls.lo.nsolves )== 3

  # test transpose solve
  x = lo.A.'\b
  x2 = zeros(x)
  linearSolveTranspose(ls, b, x2)
  @test isapprox( norm(x2 - x), 0.0) atol=1e-12
  @test ( ls.lo.nfactorizations )== 1
  @test ( ls.lo.ntsolves )== 3

  free(ls)

  return nothing
end

function test_petscmat(mesh, sbp, eqn, opts)

  println("Testing Petsc matrix-explicit linear operator")
  if PetscInitialized() == 0
    PetscInitialize()
  end
  PetscSetOptions(opts["petsc_options"])

  vals = rand(mesh.numDof, mesh.numDof)
  for i=1:mesh.numDof
    vals[i, i] = 5  # make nearly diagonally dominant
  end

  pc = PetscMatPC(mesh, sbp, eqn, opts)
  lo = PetscMatLO(pc, mesh, sbp, eqn, opts)
  ls = StandardLinearSolver(pc, lo, eqn.comm, opts)
  ctx_residual = (evalResidual,)
  t = 0.0

  # for testing
  MatSetOption(lo.A, PETSc2.MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)

  calcLinearOperator(ls, mesh, sbp, eqn, opts, ctx_residual, t)
  idx = collect(PetscInt, 1:mesh.numDof)
  idy = copy(idx)
  set_values1!(lo.A, idx, idy, vals, ADD_VALUES)

  b = rand(mesh.numDof)

  @test ( ls.shared_mat )== true
  # test product
  c = vals*b
  c2 = zeros(c)
  applyLinearOperator(ls.lo, mesh, sbp, eqn, opts, ctx_residual, t, b, c2)
  @test isapprox( norm(c - c2), 0.0) atol=1e-12

  # test transpose product
  c = vals.'*b
  applyLinearOperatorTranspose(ls.lo, mesh, sbp, eqn, opts, ctx_residual, t, b, c2)
  @test isapprox( norm(c - c2), 0.0) atol=1e-12

  setTolerances(ls, 1e-16, 1e-16, -1, -1)
  # test solve
  x = vals\b
  x2 = zeros(x)
  linearSolve(ls, b, x2)
  @test isapprox( norm(x2 - x), 0.0) atol=1e-12

  # test transpose solve
  x = vals.'\b
  x2 = zeros(x)
  linearSolveTranspose(ls, b, x2)
  @test isapprox( norm(x2 - x), 0.0) atol=1e-12

  # test product (again)
  c = vals*b
  c2 = zeros(c)
  applyLinearOperator(ls.lo, mesh, sbp, eqn, opts, ctx_residual, t, b, c2)
  @test isapprox( norm(c - c2), 0.0) atol=1e-12
  @test ( lo.nassemblies[1] )== 1

  # test transpose product (again)
  c = vals.'*b
  applyLinearOperatorTranspose(ls.lo, mesh, sbp, eqn, opts, ctx_residual, t, b, c2)
  @test isapprox( norm(c - c2), 0.0) atol=1e-12
  @test ( lo.nassemblies[1] )== 1

  # test solve (again)
  x = vals\b
  x2 = zeros(x)
  linearSolve(ls, b, x2)
  @test isapprox( norm(x2 - x), 0.0) atol=1e-12
  @test ( lo.nassemblies[1] )== 1
  @test ( lo.nsolves )== 2

  # test transpose solve (again)
  x = vals.'\b
  x2 = zeros(x)
  linearSolveTranspose(ls, b, x2)
  @test isapprox( norm(x2 - x), 0.0) atol=1e-12
  @test ( lo.nassemblies[1] )== 1
  @test ( lo.ntsolves )== 2
  @test ( pc.nsetups )== 1

  free(ls)

  return nothing
end

function test_petscmatfree(mesh, sbp, eqn, opts)

  println("testing matrix-free operator")

  vals = rand(mesh.numDof, mesh.numDof)
  for i=1:mesh.numDof
    vals[i, i] = 20  # make really diagonally dominant
  end

  # set petsc options
  petsc_opts = opts["petsc_options"]
  petsc_opts["-pc_type"] = "shell"
  PetscSetOptions(petsc_opts)

  pc = TestMatFreePC(mesh, sbp, eqn, opts)
  lo = TestMatFreeLO(pc, mesh, sbp, eqn, opts)
  ls = StandardLinearSolver(pc, lo, eqn.comm, opts)
  ctx_residual = (evalResidual, vals)
  t = 0.0

  setTolerances(ls, 1e-16, 1e-16, -1, -1)

  calcPC(ls, mesh ,sbp, eqn, opts, ctx_residual, t)
  calcLinearOperator(ls, mesh, sbp, eqn, opts, ctx_residual, t)

  b = rand(mesh.numDof)

  # test product
  c = vals*b
  c2 = zeros(c)
  applyLinearOperator(ls.lo, mesh, sbp, eqn, opts, ctx_residual, t, b, c2)
  @test isapprox( norm(c - c2), 0.0) atol=1e-12

  # test transpose product
  c = vals.'*b
  applyLinearOperatorTranspose(ls.lo, mesh, sbp, eqn, opts, ctx_residual, t, b, c2)
  @test isapprox( norm(c - c2), 0.0) atol=1e-12

  # test PC product
  x = inv(diagm(diag(vals)))*b
  x2 = zeros(x)
  applyPC(ls, mesh, sbp, eqn, opts, t, b, x2)
  @test isapprox( norm(x2 - x), 0.0) atol=1e-12

  # test PC transpose product
  x = inv(diagm(diag(vals)))*b
  x2 = zeros(x)
  applyPCTranspose(ls, mesh, sbp, eqn, opts, t, b, x2)
  @test isapprox( norm(x2 - x), 0.0) atol=1e-12

  # test solve
  x = vals\b
  x2 = zeros(x)
  linearSolve(ls, b, x2)
  @test isapprox( norm(x2 - x), 0.0) atol=1e-12

  # test transpose solve
  x = vals.'\b
  x2 = zeros(x)
  linearSolveTranspose(ls, b, x2)
  @test isapprox( norm(x2 - x), 0.0) atol=1e-12

  ### Test everything again

  # test product
  c = vals*b
  c2 = zeros(c)
  applyLinearOperator(ls.lo, mesh, sbp, eqn, opts, ctx_residual, t, b, c2)
  @test isapprox( norm(c - c2), 0.0) atol=1e-12

  # test transpose product
  c = vals.'*b
  applyLinearOperatorTranspose(ls.lo, mesh, sbp, eqn, opts, ctx_residual, t, b, c2)
  @test isapprox( norm(c - c2), 0.0) atol=1e-12

  # test PC product
  x = inv(diagm(diag(vals)))*b
  x2 = zeros(x)
  applyPC(ls, mesh, sbp, eqn, opts, t, b, x2)
  @test isapprox( norm(x2 - x), 0.0) atol=1e-12

  # test PC transpose product
  x = inv(diagm(diag(vals)))*b
  x2 = zeros(x)
  applyPCTranspose(ls, mesh, sbp, eqn, opts, t, b, x2)
  @test isapprox( norm(x2 - x), 0.0) atol=1e-12

  # test solve
  x = vals\b
  x2 = zeros(x)
  linearSolve(ls, b, x2)
  @test isapprox( norm(x2 - x), 0.0) atol=1e-12
  @test ( lo.lo_inner.nsolves )== 2

  # test transpose solve
  x = vals.'\b
  x2 = zeros(x)
  linearSolveTranspose(ls, b, x2)
  @test isapprox( norm(x2 - x), 0.0) atol=1e-12
  @test ( lo.lo_inner.ntsolves )== 2

  free(ls)

  return nothing
end

function test_petscmat_matfree(mesh, sbp, eqn, opts)
# test matrix-free linear operator with a matrix-free PC
  println("testing matrix-free operator with matrix-explicit PC")

  vals = rand(mesh.numDof, mesh.numDof)
  for i=1:mesh.numDof
    vals[i, i] = 20  # make really diagonally dominant
  end

  # set petsc options
  petsc_opts = opts["petsc_options"]
  petsc_opts["-pc_type"] = "jacobi"
  println("petsc_opts = ", petsc_opts)
  PetscSetOptions(petsc_opts)
  PetscOptionsView()

  pc = TestMatPC(mesh, sbp, eqn, opts)
  lo = TestMatFreeLO(pc, mesh, sbp, eqn, opts)
  ls = StandardLinearSolver(pc, lo, eqn.comm, opts)
  ctx_residual = (evalResidual, vals)
  t = 0.0

  setTolerances(ls, 1e-16, 1e-16, -1, -1)

  # for testing
  MatSetOption(getBasePC(pc).A, PETSc2.MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)

  calcPC(ls, mesh ,sbp, eqn, opts, ctx_residual, t)
  calcLinearOperator(ls, mesh, sbp, eqn, opts, ctx_residual, t)

  b = rand(mesh.numDof)

  # test product
  c = vals*b
  c2 = zeros(c)
  applyLinearOperator(ls.lo, mesh, sbp, eqn, opts, ctx_residual, t, b, c2)
  @test isapprox( norm(c - c2), 0.0) atol=1e-12

  # test transpose product
  c = vals.'*b
  applyLinearOperatorTranspose(ls.lo, mesh, sbp, eqn, opts, ctx_residual, t, b, c2)
  @test isapprox( norm(c - c2), 0.0) atol=1e-12

  # test PC product
  x = inv(diagm(diag(vals)))*b
  x2 = zeros(x)
  applyPC(ls, mesh, sbp, eqn, opts, t, b, x2)
  @test isapprox( norm(x2 - x), 0.0) atol=1e-12

  # test PC transpose product
  x = inv(diagm(diag(vals)))*b
  x2 = zeros(x)
  applyPCTranspose(ls, mesh, sbp, eqn, opts, t, b, x2)
  @test isapprox( norm(x2 - x), 0.0) atol=1e-12

  # test solve
  x = vals\b
  x2 = zeros(x)
  linearSolve(ls, b, x2)
  @test isapprox( norm(x2 - x), 0.0) atol=1e-12

  # test transpose solve
  x = vals.'\b
  x2 = zeros(x)
  linearSolveTranspose(ls, b, x2)
  @test isapprox( norm(x2 - x), 0.0) atol=1e-12

  free(ls)

  return nothing
end

#------------------------------------------------------------------------------
# define PC and LO
import LinearSolvers: calcPC, applyPC, applyPCTranspose, calcLinearOperator,
                      applyLinearOperator, applyLinearOperatorTranspose

# define matrix-free precondtioner for testing
# diagonal preconditioning
mutable struct TestMatFreePC <: AbstractPetscMatFreePC
  pc_inner::PetscMatFreePC
  diag::Array{Float64, 1}  # inverse of diagonal of matrix
end

function TestMatFreePC(mesh::AbstractMesh, sbp::AbstractSBP,
                       eqn::AbstractSolutionData, opts)

  pc_inner = PetscMatFreePC(mesh, sbp, eqn, opts)
  diag = Array{Float64}(mesh.numDof)

  return TestMatFreePC(pc_inner, diag)
end

function calcPC(pc::TestMatFreePC, mesh::AbstractMesh, sbp::AbstractSBP,
                eqn::AbstractSolutionData, opts::Dict, ctx_residual, t)

  f = ctx_residual[1]
  vals = ctx_residual[2]  # for testing only

  for i=1:size(vals, 1)
    pc.diag[i] = 1/vals[i, i]
  end

  setPCCtx(pc, mesh, sbp, eqn, opts, ctx_residual, t)
  pc2 = getBasePC(pc)
  pc2.nsetups += 1

  return nothing
end

function applyPC(pc::TestMatFreePC, mesh::AbstractMesh, sbp::AbstractSBP,
                eqn::AbstractSolutionData, opts::Dict, t,
                b::AbstractVector, x::AbstractVector)

  for i=1:length(pc.diag)
    x[i] = pc.diag[i]*b[i]
  end

  pc2 = getBasePC(pc)
  pc2.napplies += 1

  return nothing
end

function applyPCTranspose(pc::TestMatFreePC, mesh::AbstractMesh,
                sbp::AbstractSBP,
                eqn::AbstractSolutionData, opts::Dict, t,
                b::AbstractVector, x::AbstractVector)

  for i=1:length(pc.diag)
    x[i] = pc.diag[i]*b[i]
  end

  pc2 = getBasePC(pc)
  pc2.ntapplies += 1

  return nothing
end

#------------------------------------------------------------------------------
# define matrix-explicit PC
mutable struct TestMatPC <: AbstractPetscMatPC
  pc_inner::PetscMatPC
end

function TestMatPC(mesh::AbstractMesh, sbp::AbstractSBP,
                       eqn::AbstractSolutionData, opts)

  pc_inner = PetscMatPC(mesh, sbp, eqn, opts)

  return TestMatPC(pc_inner)
end

function calcPC(pc::TestMatPC, mesh::AbstractMesh, sbp::AbstractSBP,
                eqn::AbstractSolutionData, opts::Dict, ctx_residual, t)

  f = ctx_residual[1]
  vals = ctx_residual[2]  # for testing only

  pc2 = getBasePC(pc)
  calcPC(pc2, mesh, sbp, eqn, opts, ctx_residual, t)

  idx = collect(PetscInt, 1:mesh.numDof)
  idy = copy(idx)
  set_values1!(pc2.A, idx, idy, vals, ADD_VALUES)

  return nothing
end

# other functions do not need to be defined for matrix-explicit PC

#------------------------------------------------------------------------------
# define LO

# define fake matrix-free linear operator for testing
mutable struct TestMatFreeLO <: AbstractPetscMatFreeLO
  lo_inner::PetscMatFreeLO
  vals::Array{Float64, 2}  # the matrix
end

function TestMatFreeLO(pc::AbstractPC, mesh::AbstractMesh, sbp::AbstractSBP,
                       eqn::AbstractSolutionData, opts)

  lo_inner = PetscMatFreeLO(pc, mesh, sbp, eqn, opts)
  vals = zeros(mesh.numDof, mesh.numDof)

  return TestMatFreeLO(lo_inner, vals)
end

function calcLinearOperator(lo::TestMatFreeLO, mesh::AbstractMesh,
                sbp::AbstractSBP,
                eqn::AbstractSolutionData, opts::Dict, ctx_residual, t)

  f = ctx_residual[1]
  vals = ctx_residual[2]
  copy!(lo.vals, vals)

  setLOCtx(lo, mesh, sbp, eqn, opts, ctx_residual, t)

  return nothing
end

function applyLinearOperator(lo::TestMatFreeLO, mesh::AbstractMesh,
                sbp::AbstractSBP,
                eqn::AbstractSolutionData, opts::Dict, ctx_residual, t,
                x::AbstractVector, b::AbstractVector)

  smallmatvec!(lo.vals, x, b)

  return nothing
end

function applyLinearOperatorTranspose(lo::TestMatFreeLO, mesh::AbstractMesh,
                sbp::AbstractSBP,
                eqn::AbstractSolutionData, opts::Dict, ctx_residual, t,
                x::AbstractVector, b::AbstractVector)

  smallmatTvec!(lo.vals, x, b)

  return nothing
end

add_func1!(AdvectionTests, test_linearsolver, [TAG_SHORTTEST])


text:
"""
    ReshapedDistribution(d::Distribution{<:ArrayLikeVariate}, dims::Dims{N})

Distribution of the `N`-dimensional random variable `reshape(X, dims)` where `X` is a random
variable wth distribution `d`.

It is recommended to not use `reshape` instead of the constructor of `ReshapedDistribution`
directly since `reshape` can return more optimized distributions for specific types of `d`
and number of dimensions `N`.
"""
struct ReshapedDistribution{N,S<:ValueSupport,D<:Distribution{<:ArrayLikeVariate,S}} <: Distribution{ArrayLikeVariate{N},S}
    dist::D
    dims::Dims{N}

    function ReshapedDistribution(dist::Distribution{<:ArrayLikeVariate,S}, dims::Dims{N}) where {N,S<:ValueSupport}
        _reshape_check_dims(dist, dims)
        return new{N,S,typeof(dist)}(dist, dims)
    end
end

function _reshape_check_dims(dist::Distribution{<:ArrayLikeVariate}, dims::Dims)
    (all(d > 0 for d in dims) && length(dist) == prod(dims)) ||
        throw(ArgumentError("dimensions $(dims) do not match size of source distribution $(size(dist))"))
end

Base.size(d::ReshapedDistribution) = d.dims
Base.eltype(::Type{ReshapedDistribution{<:Any,<:ValueSupport,D}}) where {D} = eltype(D)

partype(d::ReshapedDistribution) = partype(d.dist)
params(d::ReshapedDistribution) = (d.dist, d.dims)

function insupport(d::ReshapedDistribution{N}, x::AbstractArray{<:Real,N}) where {N}
    return size(d) == size(x) && insupport(d.dist, reshape(x, size(d.dist)))
end

mean(d::ReshapedDistribution) = reshape(mean(d.dist), size(d))
var(d::ReshapedDistribution) = reshape(var(d.dist), size(d))
cov(d::ReshapedDistribution) = reshape(cov(d.dist), length(d), length(d))
function cov(d::ReshapedDistribution{2}, ::Val{false})
    n, p = size(d)
    return reshape(cov(d), n, p, n, p)
end

mode(d::ReshapedDistribution) = reshape(mode(d.dist), size(d))

# TODO: remove?
rank(d::ReshapedDistribution{2}) = minimum(size(d))

# logpdf evaluation
# have to fix method ambiguity due to default fallback for `MatrixDistribution`...
_logpdf(d::ReshapedDistribution{N}, x::AbstractArray{<:Real,N}) where {N} = __logpdf(d, x)
_logpdf(d::ReshapedDistribution{2}, x::AbstractMatrix{<:Real}) = __logpdf(d, x)
function __logpdf(d::ReshapedDistribution{N}, x::AbstractArray{<:Real,N}) where {N}
    dist = d.dist
    return @inbounds logpdf(dist, reshape(x, size(dist)))
end

# loglikelihood
# useful if the original distribution defined more optimized methods
@inline function loglikelihood(
    d::ReshapedDistribution{N},
    x::AbstractArray{<:Real,N},
) where {N}
    @boundscheck begin
        size(x) == size(d) ||
            throw(DimensionMismatch("inconsistent array dimensions"))
    end
    dist = d.dist
    return @inbounds loglikelihood(dist, reshape(x, size(dist)))
end
@inline function loglikelihood(
    d::ReshapedDistribution{N},
    x::AbstractArray{<:Real,M},
) where {N,M}
    @boundscheck begin
        M > N ||
            throw(DimensionMismatch(
                "number of dimensions of `x` ($M) must be greater than number of dimensions of `d` ($N)"
            ))
        ntuple(i -> size(x, i), Val(N)) == size(d) ||
            throw(DimensionMismatch("inconsistent array dimensions"))
    end
    dist = d.dist
    trailingsize = ntuple(i -> size(x, N + i), Val(M - N))
    return @inbounds loglikelihood(dist, reshape(x, size(dist)..., trailingsize...))
end

# sampling
function _rand!(
    rng::AbstractRNG,
    d::ReshapedDistribution{N},
    x::AbstractArray{<:Real,N}
) where {N}
    dist = d.dist
    @inbounds rand!(rng, dist, reshape(x, size(dist)))
    return x
end

"""
    reshape(d::Distribution{<:ArrayLikeVariate}, dims::Int...)
    reshape(d::Distribution{<:ArrayLikeVariate}, dims::Dims)

Return a [`Distribution`](@ref) of `reshape(X, dims)` where `X` is a random variable with
distribution `d`.

The default implementation returns a [`ReshapedDistribution`](@ref). However, it can return
more optimized distributions for specific types of distributions and numbers of dimensions.
Therefore it is recommended to use `reshape` instead of the constructor of
`ReshapedDistribution`.

# Implementation

Since `reshape(d, dims::Int...)` calls `reshape(d, dims::Dims)`, one should implement
`reshape(d, ::Dims)` for desired distributions `d`.

See also: [`vec`](@ref)
"""
function Base.reshape(dist::Distribution{<:ArrayLikeVariate}, dims::Dims)
    return ReshapedDistribution(dist, dims)
end
function Base.reshape(dist::Distribution{<:ArrayLikeVariate}, dims1::Int, dims::Int...)
    return reshape(dist, (dims1, dims...))
end

"""
    vec(d::Distribution{<:ArrayLikeVariate})

Return a [`MultivariateDistribution`](@ref) of `vec(X)` where `X` is a random variable with
distribution `d`.

The default implementation returns a [`ReshapedDistribution`](@ref). However, it can return
more optimized distributions for specific types of distributions and numbers of dimensions.
Therefore it is recommended to use `vec` instead of the constructor of
`ReshapedDistribution`.

# Implementation

Since `vec(d)` is defined as `reshape(d, length(d))` one should implement
`reshape(d, ::Tuple{Int})` rather than `vec`.

See also: [`reshape`](@ref)
"""
Base.vec(dist::Distribution{<:ArrayLikeVariate}) = reshape(dist, length(dist))

# avoid unnecessary wrappers
function Base.reshape(
    dist::ReshapedDistribution{<:Any,<:ValueSupport,<:MultivariateDistribution},
    dims::Tuple{Int},
)
    _reshape_check_dims(dist, dims)
    return dist.dist
end

function Base.reshape(dist::MultivariateDistribution, dims::Tuple{Int})
    _reshape_check_dims(dist, dims)
    return dist
end

# specialization for flattened `MatrixNormal`
function Base.reshape(dist::MatrixNormal, dims::Tuple{Int})
    _reshape_check_dims(dist, dims)
    return MvNormal(vec(dist.M), kron(dist.V, dist.U))
end


text:
##
# ProductFun represents f(x,y) by Fun(.coefficients[k](x),.space[2])(y)
# where all coefficients are in the same space
##

export ProductFun

struct ProductFun{S<:UnivariateSpace,V<:UnivariateSpace,SS<:AbstractProductSpace,T} <: BivariateFun{T}
    coefficients::Vector{VFun{S,T}}     # coefficients are in x
    space::SS
end

ProductFun(cfs::Vector{VFun{S,T}},sp::AbstractProductSpace{Tuple{S,V},DD}) where {S<:UnivariateSpace,V<:UnivariateSpace,T<:Number,DD} =
    ProductFun{S,V,typeof(sp),T}(cfs,sp)
ProductFun(cfs::Vector{VFun{S,T}},sp::AbstractProductSpace{Tuple{W,V},DD}) where {S<:UnivariateSpace,V<:UnivariateSpace,
         W<:UnivariateSpace,T<:Number,DD} =
   ProductFun{W,V,typeof(sp),T}(VFun{W,T}[Fun(cfs[k],columnspace(sp,k)) for k=1:length(cfs)],sp)

Base.size(f::ProductFun,k::Integer) =
    k==1 ? mapreduce(ncoefficients,max,f.coefficients) : length(f.coefficients)
Base.size(f::ProductFun) = (size(f,1),size(f,2))

## Construction in an AbstractProductSpace via a Matrix of coefficients

function ProductFun(cfs::AbstractMatrix{T},sp::AbstractProductSpace{Tuple{S,V},DD};
  tol::Real=100eps(T),chopping::Bool=false) where {S<:UnivariateSpace,V<:UnivariateSpace,T<:Number,DD}
    if chopping
        ncfs,kend=norm(cfs,Inf),size(cfs,2)
        if kend > 1 while isempty(chop(cfs[:,kend],ncfs*tol)) kend-=1 end end
        ret=VFun{S,T}[Fun(columnspace(sp,k),chop(cfs[:,k],ncfs*tol)) for k=1:max(kend,1)]
        ProductFun{S,V,typeof(sp),T}(ret,sp)
    else
        ret=VFun{S,T}[Fun(columnspace(sp,k),cfs[:,k]) for k=1:size(cfs,2)]
        ProductFun{S,V,typeof(sp),T}(ret,sp)
    end
end

## Construction in a ProductSpace via a Vector of Funs

function ProductFun(M::Vector{VFun{S,T}},dy::V) where {S<:UnivariateSpace,V<:UnivariateSpace,T<:Number}
    funs=VFun{S,T}[Mk for Mk in M]
    ProductFun{S,V,ProductSpace{S,V},T}(funs,ProductSpace(S[space(fun) for fun in funs],dy))
end

## Adaptive construction

function ProductFun(f::Function,sp::AbstractProductSpace{Tuple{S,V}};tol=100eps()) where {S<:UnivariateSpace,V<:UnivariateSpace}
    for n = 50:100:5000
        X = coefficients(ProductFun(dynamic(f),sp,n,n;tol=tol))
        if size(X,1)<n && size(X,2)<n
            return ProductFun(X,sp;tol=tol)
        end
    end
    @warn "Maximum grid size of ("*string(5000)*","*string(5000)*") reached"
    ProductFun(f,sp,5000,5000;tol=tol,chopping=true)
end

## ProductFun values to coefficients

function ProductFun(f::Function,S::AbstractProductSpace,M::Integer,N::Integer;tol=100eps())
    xy = checkpoints(S)
    T = promote_type(eltype(f(first(xy)...)),rangetype(S))
    ptsx,ptsy=points(S,M,N)
    vals=T[f(ptsx[k,j],ptsy[k,j]) for k=1:size(ptsx,1), j=1:size(ptsx,2)]
    ProductFun(transform!(S,vals),S;tol=tol,chopping=true)
end
ProductFun(f::Function,S::TensorSpace) = ProductFun(LowRankFun(dynamic(f),S))

ProductFun(f,dx::Space,dy::Space)=ProductFun(dynamic(f),TensorSpace(dx,dy))

ProductFun(f::Function,dx::Space,dy::Space)=ProductFun(dynamic(f),TensorSpace(dx,dy))

## Domains promoted to Spaces

ProductFun(f::Function,D::Domain,M::Integer,N::Integer) = ProductFun(dynamic(f),Space(D),M,N)
ProductFun(f::Function,d::Domain) = ProductFun(dynamic(f),Space(d))
ProductFun(f::Function, dx::Domain{<:Number}, dy::Domain{<:Number}) = ProductFun(dynamic(f),Space(dx),Space(dy))
ProductFun(f::Function) = ProductFun(dynamic(f),ChebyshevInterval(),ChebyshevInterval())

## Conversion from other 2D Funs

ProductFun(f::LowRankFun) = ProductFun(coefficients(f),space(f,1),space(f,2))
ProductFun(f::Fun{S}) where {S<:AbstractProductSpace} = ProductFun(coefficientmatrix(f),space(f))

## Conversion to other ProductSpaces with the same coefficients

ProductFun(f::ProductFun,sp::TensorSpace)=space(f)==sp ? f : ProductFun(coefficients(f,sp),sp)
ProductFun(f::ProductFun{S,V,SS},sp::ProductDomain) where {S,V,SS<:TensorSpace}=ProductFun(f,Space(sp))

function ProductFun(f::ProductFun,sp::AbstractProductSpace)
    u=Array{VFun{typeof(columnspace(sp,1)),cfstype(f)}}(length(f.coefficients))

    for k=1:length(f.coefficients)
        u[k]=Fun(f.coefficients[k],columnspace(sp,k))
    end

    ProductFun(u,sp)
end

## For specifying spaces by anonymous function

ProductFun(f::Function,SF::Function,T::Space,M::Integer,N::Integer) =
    ProductFun(dynamic(f),typeof(SF(1))[SF(k) for k=1:N],T,M)

## Conversion of a constant to a ProductFun

ProductFun(c::Number,sp::BivariateSpace) = ProductFun([Fun(c,columnspace(sp,1))],sp)
ProductFun(f::Fun,sp::BivariateSpace) = ProductFun([Fun(f,columnspace(sp,1))],sp)

## Utilities

function funlist2coefficients(f::Vector{VFun{S,T}}) where {S,T}
    A=zeros(T,mapreduce(ncoefficients,max,f),length(f))
    for k=1:length(f)
        A[1:ncoefficients(f[k]),k]=f[k].coefficients
    end
    A
end

function pad(f::ProductFun{S,V,SS,T},n::Integer,m::Integer) where {S,V,SS,T}
    ret=Array{VFun{S,T}}(undef, m)
    cm=min(length(f.coefficients),m)
    for k=1:cm
        ret[k]=pad(f.coefficients[k],n)
    end

    for k=cm+1:m
        ret[k] = zeros(columnspace(f,k))
    end
    ProductFun{S,V,SS,T}(ret,f.space)
end

function pad!(f::ProductFun{S,V,SS,T},::Colon,m::Integer) where {S,V,SS,T}
    cm=length(f.coefficients)
    resize!(f.coefficients,m)

    for k=cm+1:m
        f.coefficients[k]=zeros(columnspace(f,k))
    end
    f
end

coefficients(f::ProductFun)=funlist2coefficients(f.coefficients)

function coefficients(f::ProductFun,ox::Space,oy::Space)
    T=cfstype(f)
    m=size(f,1)
    B=Matrix{T}(m,length(f.coefficients))
    # convert in x direction
    #TODO: adaptively grow in x?
    for k=1:length(f.coefficients)
        B[:,k]=pad!(coefficients(f.coefficients[k],ox),m)
    end

    # convert in y direction
    for k=1:size(B,1)
        ccfs=coefficients(vec(B[k,:]),space(f,2),oy)
        if length(ccfs)>size(B,2)
            B=pad(B,size(B,1),length(ccfs))
        end
        B[k,1:length(ccfs)]=ccfs
        for j=length(ccfs)+1:size(B,2)
            B[k,j]=zero(T)
        end
    end

    B
end

(f::ProductFun)(x,y) = evaluate(f,x,y)
(f::ProductFun)(x,y,z) = evaluate(f,x,y,z)

coefficients(f::ProductFun,ox::TensorSpace) = coefficients(f,ox[1],ox[2])

values(f::ProductFun{S,V,SS,T}) where {S,V,SS,T} = itransform!(space(f),coefficients(f))

vecpoints(f::ProductFun{S,V,SS},k) where {S,V,SS<:TensorSpace} = points(f.space[k],size(f,k))

space(f::ProductFun) = f.space
columnspace(f::ProductFun,k) = columnspace(space(f),k)

domain(f::ProductFun) = domain(f.space)
#domain(f::ProductFun,k)=domain(f.space,k)
canonicaldomain(f::ProductFun) = canonicaldomain(space(f))

function canonicalevaluate(f::ProductFun{S,V,SS,T},x::Number,::Colon) where {S,V,SS,T}
    cd = canonicaldomain(f)
    Fun(setdomain(factor(space(f),2),factor(cd,2)),
                    T[setdomain(fc,factor(cd,1))(x) for fc in f.coefficients])
end
canonicalevaluate(f::ProductFun,x::Number,y::Number) = canonicalevaluate(f,x,:)(y)
canonicalevaluate(f::ProductFun{S,V,SS},x::Colon,y::Number) where {S,V,SS<:TensorSpace} =
    evaluate(transpose(f),y,:)  # doesn't make sense For general product fon without specifying space

canonicalevaluate(f::ProductFun,xx::AbstractVector,yy::AbstractVector) =
    transpose(hcat([evaluate(f,x,:)(yy) for x in xx]...))

evaluate(f::ProductFun,x,y) = canonicalevaluate(f,tocanonical(f,x,y)...)
evaluate(f::ProductFun,x,y,z) = canonicalevaluate(f,tocanonical(f,x,y,z)...)

# TensorSpace does not use map
evaluate(f::ProductFun{S,V,SS,T},x::Number,::Colon) where {S<:UnivariateSpace,V<:UnivariateSpace,SS<:TensorSpace,T} =
    Fun(factor(space(f),2),T[g(x) for g in f.coefficients])

evaluate(f::ProductFun{S,V,SS,T},x::Number,y::Number) where {S<:UnivariateSpace,V<:UnivariateSpace,SS<:TensorSpace,T} =
    evaluate(f,x,:)(y)

evaluate(f::ProductFun,x) = evaluate(f,x...)

*(c::Number,f::F) where {F<:ProductFun} = F(c*f.coefficients,f.space)
*(f::ProductFun,c::Number) = c*f

function chop(f::ProductFun{S},es...) where S
    kend=size(f,2)
    while kend > 1 && isempty(chop(f.coefficients[kend].coefficients,es...))
        kend-=1
    end
    ret=VFun{S,cfstype(f)}[Fun(space(f.coefficients[k]),chop(f.coefficients[k].coefficients,es...)) for k=1:max(kend,1)]

    typeof(f)(ret,f.space)
end

##TODO: following assumes f is never changed....maybe should be deepcopy?
function +(f::F,c::Number) where F<:ProductFun
    cfs=copy(f.coefficients)
    cfs[1]+=c
    F(cfs,f.space)
end
+(c::Number,f::ProductFun) = f+c
-(f::ProductFun,c::Number) = f+(-c)
-(c::Number,f::ProductFun) = c+(-f)

function +(f::ProductFun,g::ProductFun)
    if f.space == g.space
        if size(f,2) >= size(g,2)
            @assert f.space==g.space
            cfs = copy(f.coefficients)
            for k=1:size(g,2)
                cfs[k]+=g.coefficients[k]
            end

            ProductFun(cfs,f.space)
        else
            g+f
        end
    else
        s=conversion_type(f.space,g.space)
        ProductFun(f,s)+ProductFun(g,s)
    end
end

-(f::ProductFun) = (-1)*f
-(f::ProductFun,g::ProductFun) = f+(-g)

*(B::Fun,f::ProductFun) = ProductFun(map(c->B*c,f.coefficients),space(f))
*(f::ProductFun,B::Fun) = transpose(B*transpose(f))

LowRankFun(f::ProductFun{S,V,SS}) where {S,V,SS<:TensorSpace} = LowRankFun(f.coefficients,factor(space(f),2))
LowRankFun(f::Fun) = LowRankFun(ProductFun(f))

function differentiate(f::ProductFun{S,V,SS},j::Integer) where {S,V,SS<:TensorSpace}
    if j==1
        df=map(differentiate,f.coefficients)
        ProductFun(df,space(first(df)),factor(space(f),2))
    else
        transpose(differentiate(transpose(f),1))
    end
end

# If the transpose of the space exists, then the transpose of the ProductFun exists
Base.transpose(f::ProductFun{S,V,SS,T}) where {S,V,SS,T} =
    ProductFun(transpose(coefficients(f)),transpose(space(f)))

for op in (:(Base.sin),:(Base.cos))
    @eval ($op)(f::ProductFun) =
        Fun(space(f),transform!(space(f),$op(values(pad(f,size(f,1)+20,size(f,2))))))
end

^(f::ProductFun,k::Integer) =
    Fun(space(f),transform!(space(f),values(pad(f,size(f,1)+20,size(f,2))).^k))

for op = (:(Base.real),:(Base.imag),:(Base.conj))
    @eval ($op)(f::ProductFun{S,V,SS}) where {S,V<:RealSpace,SS<:TensorSpace} =
        ProductFun(map($op,f.coefficients),space(f))
end

#For complex bases
Base.real(f::ProductFun{S,V,SS}) where {S,V,SS<:TensorSpace} =
    transpose(real(transpose(ProductFun(real(u.coefficients),space(u)))))-transpose(imag(transpose(ProductFun(imag(u.coefficients),space(u)))))
Base.imag(f::ProductFun{S,V,SS}) where {S,V,SS<:TensorSpace} =
    transpose(real(transpose(ProductFun(imag(u.coefficients),space(u)))))+transpose(imag(transpose(ProductFun(real(u.coefficients),space(u)))))

## Call LowRankFun version
# TODO: should cumsum and integrate return TensorFun or lowrankfun?
for op in (:(Base.sum),:(Base.cumsum),:integrate)
    @eval $op(f::ProductFun{S,V,SS},n...) where {S,V,SS<:TensorSpace} = $op(LowRankFun(f),n...)
end

## ProductFun transform

# function transform{ST<:Space,N<:Number}(::Type{N},S::Vector{ST},T::Space,V::AbstractMatrix)
#     @assert length(S)==size(V,2)
#     # We assume all S spaces have same domain/points
#     C=Vector{N}(size(V)...)
#     for k=1:size(V,1)
#         C[k,:]=transform(T,vec(V[k,:]))
#     end
#     for k=1:size(C,2)
#         C[:,k]=transform(S[k],C[:,k])
#     end
#     C
# end
# transform{ST<:Space,N<:Real}(S::Vector{ST},T::Space{Float64},V::AbstractMatrix{N})=transform(Float64,S,T,V)
# transform{ST<:Space}(S::Vector{ST},T::Space,V::AbstractMatrix)=transform(Complex{Float64},S,T,V)

for op in (:tocanonical,:fromcanonical)
    @eval $op(f::ProductFun,x...) = $op(space(f),x...)
end


text:
using Test
using NumericalTricks
using StaticArrays
using LinearAlgebra
using Rotations
using ForwardDiff: Partials, Dual, value

include("test_exports.jl")
include("test_div_by_zero.jl")
# include("test_matrix_factor_derivatives.jl")
include("test_basic_dh.jl")
include("test_geometry_kernel.jl")
include("test_poly_approx.jl")
include("test_utility.jl")
include("test_vector_projections.jl")
include("test_rotations.jl")


text:
module SIMDPirates

using MacroTools: postwalk, @capture
using VectorizationBase
using VectorizationBase: llvmtype, AbstractSIMDVector, AbstractStructVec, AbstractStridedPointer

export  Vec, SVec, VE,
    @pirate,
    extract_data,
    vbroadcast,
    vload,
    vloada,
    vstore!,
    vstorea!,
    shufflevector,
    vifelse,
    vfma, vmuladd,
    vsqrt, rsqrt, vinv, vabs2,
    vadd, vsub, vmul, vfdiv, vsum, vprod,
    vfmadd, vfnmadd, vfmsub, vfnmsub,
    gather, scatter!

function vector_args(args)
    vecargs =       [isa(arg, Symbol) ? :($arg::Vec{N,T})             : arg for arg ∈ args]
    abstractargs =  [isa(arg, Symbol) ? :($arg::AbstractSIMDVector{N,T})    : arg for arg ∈ args]
    structargs =    [isa(arg, Symbol) ? :($arg::AbstractStructVec{N,T})     : arg for arg ∈ args]
    vecargs, abstractargs, structargs
end

function parse_func_decl(expr::Expr)
    @assert (expr.head === :(=) || expr.head === :function) "Failed to match function declaration: $expr."
    w = expr.args[1]
    body = expr.args[2]
    WT = @view(w.args[2:end])
    func = first(w.args)
    f = first(func.args)
    args = @view(func.args[2:end])
    f, args, WT, body
end
macro vectordef(rename, expr)
    f, args, R, body = parse_func_decl(expr)
    vecargs, abstractargs, structargs = vector_args(args)
    q = quote
        @inline $rename($(vecargs...)) where {$(R...)} = $body
        @inline $rename($(abstractargs...)) where {$(R...)} = SVec($body)
        @inline $f($(structargs...)) where {$(R...)} = SVec($body)
    end
    esc(q)
end
macro evectordef(rename, expr)
    f, args, R, body = parse_func_decl(expr)
    vecargs, abstractargs, structargs = vector_args(args)
    q = quote
        @inline $rename($(vecargs...)) where {$(R...)} = $body
        @inline $rename($(abstractargs...)) where {$(R...)} = SVec($body)
    end
    esc(q)
end

include("type_definitions.jl")
include("llvm_utils.jl")
include("llvmwrap.jl")
include("conditionals.jl")
include("integer_arithmetic.jl")
include("floating_point_arithmetic.jl")
include("memory.jl")
include("shufflevector.jl")
include("special.jl")
include("pirate.jl")
include("precompile.jl")
_precompile_()

function __init__()
    _precompile_()
end

end # module


text:
#=
 Jaywalking model in WebPPL.

  This model simulates how many people (and why) that will jaywalk
  of a crossing line with red light, given that they have some
  prospensity of jaywalk if some <limit> number of people already
  has jaywalked .

  The mode:
     1) Assign some number n (might be fixed)
     2) Foreach agent: generate a random tendency to jaywalk.
     3) For each time slot t=0..n-1
          foreach a=agent (0..n-1):
            if there are more than or equal number of <limit> number
            of people that has already jaywalked, then a will also
            jaywalk.

  The number of people that - eventually - jaywalks depends very
  much on the distribution of the limits and espespecially
  the probability of the 0 limits jaywalkers, i.e. those that always
  jaywalk.

  Here is an some example where n_max = 15 (
 - endTime: number of timeslots
 - n=number of people
 - t0: number of people that initially jaywalk
 - tn: number of people that eventually will jaywalk

  Mean n:  8.5141
  Distributions of variable n (num:0)
  3.00000 =>     762  (0.076200)
  13.00000 =>     740  (0.074000)
  15.00000 =>     728  (0.072800)
  8.00000 =>     721  (0.072100)
  5.00000 =>     720  (0.072000)
  4.00000 =>     717  (0.071700)
  9.00000 =>     715  (0.071500)
  11.00000 =>     714  (0.071400)
  14.00000 =>     714  (0.071400)
  12.00000 =>     712  (0.071200)
  6.00000 =>     702  (0.070200)
  10.00000 =>     696  (0.069600)
  2.00000 =>     687  (0.068700)
  7.00000 =>     672  (0.067200)

  Mean t0:  0.5049
  Distributions of variable t0 (num:0)
  0.00000 =>    5896  (0.589600)
  1.00000 =>    3272  (0.327200)
  2.00000 =>     729  (0.072900)
  3.00000 =>      95  (0.009500)
  4.00000 =>       6  (0.000600)
  5.00000 =>       2  (0.000200)

  Mean tn:  6.3268
  Distributions of variable tn (num:0)
  0.00000 =>    5896  (0.589600)
  15.00000 =>     661  (0.066100)
  16.00000 =>     577  (0.057700)
  17.00000 =>     512  (0.051200)
  14.00000 =>     491  (0.049100)
  13.00000 =>     403  (0.040300)
  18.00000 =>     393  (0.039300)
  19.00000 =>     257  (0.025700)
  12.00000 =>     254  (0.025400)
  11.00000 =>     176  (0.017600)
  20.00000 =>     133  (0.013300)
  10.00000 =>      74  (0.007400)
  21.00000 =>      64  (0.006400)
  22.00000 =>      34  (0.003400)
  9.00000 =>      28  (0.002800)
  8.00000 =>      21  (0.002100)
  23.00000 =>      15  (0.001500)
  7.00000 =>       5  (0.000500)
  6.00000 =>       3  (0.000300)
  24.00000 =>       2  (0.000200)
  5.00000 =>       1  (0.000100)

  Cf ~/webppl/jaywalking.wppl

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function jaywalking(n_max=15)
    endTime = n_max*2

    # n = 15; # number of people
    n ~ DiscreteUniform(2,n_max) # number of people

    # The limits for jaywalking, i.e.
    # how many persons that have to jaywalk
    # before this agent jaywlk.
    limit ~ filldist(DiscreteUniform(0,n*2-1),endTime)

    t0 ~ Dirac(sum([limit[i] == 0 ? 1 : 0 for i in 1:n]))

    # The number of people that will jaywalk at time i
    t = tzeros(endTime)
    for i in 1:endTime
        t[i] ~ Dirac(sum([ t0 > 0 && limit[j] <= i ? 1 : 0 for j in 1:endTime]))
    end
    tn ~ Dirac(t0 > 0  ? t[n-1] : 0)
end

model = jaywalking()

chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 100_000)
# chns = sample(model, PG(5), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

# display(chns)
# display(plot(chns))

println("Mean n:  ", mean_val(chns, :n))
show_var_dist_pct(chns, :n)
println("Mean t0:  ", mean_val(chns, :t0))
show_var_dist_pct(chns, :t0)
println("Mean tn:  ", mean_val(chns, :tn))
show_var_dist_pct(chns, :tn)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
## LICENSE
##   Copyright (c) 2018 GEOPHYSTECH LLC
##
##   Licensed under the Apache License, Version 2.0 (the "License");
##   you may not use this file except in compliance with the License.
##   You may obtain a copy of the License at
##
##       http://www.apache.org/licenses/LICENSE-2.0
##
##   Unless required by applicable law or agreed to in writing, software
##   distributed under the License is distributed on an "AS IS" BASIS,
##   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
##   See the License for the specific language governing permissions and
##   limitations under the License.

## initial release by Andrey Stepnov a.stepnov@geophsytech.ru

"""
**ON GRID**

`gmpe_as2008(eq::Earthquake,config_as2008::Params_as2008,grid::Array{Point_vs30};min_val::Number)`

where `min_val=0` by default

Output will be 1-d `Array{Point_pga_out}` with points based on input grid with `pga > min_val` (pga, is Acceleration of gravity in percent (%g) rounded to ggg.gg)

**without grid**

`gmpe_as2008(eq::Earthquake,config::Params_as2008,VS30::Number=350,distance::Int64=1000)`

where `VS30=30` [m/s^2], `distance=1000` [km] by default.

Output will be 1-d `Array{Float64}` with `1:distance` values of `pga` (that is Acceleration of gravity in percent (%g) rounded to ggg.gg)

**EXAMPLES:**
```
gmpe_as2008(eq,config_as2008,grid) # for PGA on GRID
gmpe_as2008(eq,config_as2008) # for without input grid
```

**Model parameters**

Please, see `examples/as-2008.conf`
"""
## AS2008 PGA modeling ON GRID
function gmpe_as2008(eq::Earthquake,config::Params_as2008,grid::Array{Point_vs30};min_val::Number=0)
  vs30_row_num = length(grid[:,1])
  eq.moment_mag == 0 ? magnitude = eq.local_mag : magnitude = eq.moment_mag
  epicenter = LatLon(eq.lat, eq.lon)
  # define t6
  if magnitude < 5.5
    t6 = 1
  elseif magnitude <= 6.5 && magnitude >= 5.5
    t6 = (0.5 * (6.5 - magnitude) + 0.5)
  else
    t6 = 0.5
  end
  # modeling
  if config.ground_motion_type == "PGA"
    output_data = Array{Point_pga_out}(undef,0)
    out_type = Point_pga_out
  end
  for i=1:vs30_row_num
    # rrup
    current_point = LatLon(grid[i].lat,grid[i].lon)
    r_rup = sqrt((distance(current_point,epicenter)/1000)^2 + eq.depth^2)
    # F1
    if magnitude <= config.c1
      f1 = config.a1 + config.a4 * (magnitude - config.c1) +
        config.a8 * (8.5 - magnitude)^2 +
        (config.a2 + config.a3 * (magnitude - config.c1)) *
        log(sqrt(r_rup^2 + config.c4^2))
    else
      f1 = config.a1 + config.a5 * (magnitude - config.c1) +
        config.a8 * (8.5 - magnitude)^2 +
        (config.a2 + config.a3 * (magnitude - config.c1)) *
            log(sqrt(r_rup^2 + config.c4^2))
    end
    # F8
    if r_rup < 100
      f8 = 0;
    else
      f8 = config.a18 * (r_rup - 100) * t6
    end
    # PGA1100
    pga1100 = exp((f1 + (config.a10 + config.b * config.n) *
               log(config.vs30_1100 / config.vlin) + f8))
    # F5
    if grid[i].vs30 < config.vlin
      f5 =  config.a10 * log(grid[i].vs30 / config.vlin) -
        config.b * log(pga1100 + config.c) + config.b *
        log(pga1100 + config.c * (grid[i].vs30 / config.vlin)^config.n)
    elseif (grid[i].vs30 > config.vlin) && (grid[i].vs30 < config.v1)
      f5 = (config.a10 + config.b * config.n) *
        log(grid[i].vs30 / config.vlin)
    else
      f5 = (config.a10 + config.b * config.n) *
        log(config.v1 / config.vlin)
    end
    if config.ground_motion_type == "PGA"
      motion = round((exp(f1 + f5 + f8) * 100),digits=2)
    end
    if motion >= min_val
      output_data = push!(output_data, out_type(grid[i].lon,grid[i].lat,motion))
    end
  end
  return output_data
end
## AS2008 PGA modeling for PLOTTING
function gmpe_as2008(eq::Earthquake,config::Params_as2008;VS30::Number=350,distance::Number=1000)
  eq.moment_mag == 0 ? magnitude = eq.local_mag : magnitude = eq.moment_mag
  # define t6
  if magnitude < 5.5
    t6 = 1
  elseif magnitude <= 6.5 && magnitude >= 5.5
    t6 = (0.5 * (6.5 - magnitude) + 0.5)
  else
    t6 = 0.5
  end
  output_data = Array{Float64}(undef,0)
  for i=1:distance
    # rrup
    r_rup = sqrt((i)^2 + eq.depth^2)
    # F1
    if magnitude <= config.c1
      f1 = config.a1 + config.a4 * (magnitude - config.c1) +
        config.a8 * (8.5 - magnitude)^2 +
        (config.a2 + config.a3 * (magnitude - config.c1)) *
        log(sqrt(r_rup^2 + config.c4^2))
    else
      f1 = config.a1 + config.a5 * (magnitude - config.c1) +
        config.a8 * (8.5 - magnitude)^2 +
        (config.a2 + config.a3 * (magnitude - config.c1)) *
            log(sqrt(r_rup^2 + config.c4^2))
    end
    # F8
    if r_rup < 100
      f8 = 0;
    else
      f8 = config.a18 * (r_rup - 100) * t6
    end
    # PGA1100
    pga1100 = exp((f1 + (config.a10 + config.b * config.n) *
               log(config.vs30_1100 / config.vlin) + f8))
    # F5
    if VS30 < config.vlin
      f5 =  config.a10 * log(VS30 / config.vlin) -
        config.b * log(pga1100 + config.c) + config.b *
        log(pga1100 + config.c * (VS30 / config.vlin)^config.n)
    elseif (VS30 > config.vlin) && (VS30 < config.v1)
      f5 = (config.a10 + config.b * config.n) *
        log(VS30 / config.vlin)
    else
      f5 = (config.a10 + config.b * config.n) *
        log(config.v1 / config.vlin)
    end
    if config.ground_motion_type == "PGA"
      motion = round((exp(f1 + f5 + f8) * 100),digits=2)
    end
    output_data = push!(output_data, motion)
  end
  return output_data
end


text:
export DynamicAlgebraic

struct DynamicAlgebraic
    dynamic::Function
    algebraic::Function
end


text:
using Revise
import MitosisStochasticDiffEq as MSDE
using Test, Random
using LinearAlgebra
using BenchmarkTools
using MitosisStochasticDiffEq: NoiseGrid, EM

#=
function MSDE.saveit!(uu::Vector, u::Tuple{Int64, Float64, Vector{Float64}}, P)
  push!(uu, (u[1], u[2], copy(u[3])))
end
=#

# define SDE function
f(u,p,t) = p[1]*u + p[2]
g(u,p,t) = p[3]*u

function f!(du, u, p, t)
  mul!(du, p[1], u)
  @. du += p[2]
  du
end

g!(du, u, p, t) = (@. du = p[3]*u)
function gstep!(dx, _, u, p, t, dw, _)
  @. dx += (p[3]*u)*dw
end
# time span
tstart = 0.0
tend = 1.0
dt = 0.00001
trange = tstart:dt:tend

B, β, σ̃ = -fill(0.1,1,1), [0.2], 1.3
p = (B, β, σ̃)

u0 = rand(1)

k_oop = MSDE.SDEKernel(f, g, trange, p)
k_iip = MSDE.SDEKernel(f!, g!, trange, p)
k_iip2 = MSDE.SDEKernel!(f!, g!, gstep!, trange, p; ws = copy(u0))

# P manually
struct customP{θType}
  θ::θType
end

function MSDE.dZ!(u, dz, Z, P::customP)
  i = u[1]
  dw = dz[3]
  @. dw = Z[i+1][3] -  Z[i][3]
  (Z[i+1][1] - Z[i][1], Z[i+1][2] - Z[i][2], dw)
end

# mul!(C, A, B, α, β) -> C = A B α + C β.

@inline function MSDE.tangent!(du, u, dz, P::customP)
  x = du[3]
  mul!(x, P.θ[1], u[3], dz[2], false)
  @. x += P.θ[2]*dz[2] + P.θ[3]*u[3]*dz[3]
  (dz[1], dz[2], x)
end

@inline function MSDE.exponential_map!(u, du, P::customP)
  x = u[3]
  @. x += du[3]
  (u[1] + du[1], u[2] + du[2], x)
end

# pass noise process and compare with EM()
Ws = cumsum([[zero(u0)];[sqrt(trange[i+1]-ti)*randn(size(u0))
        for (i,ti) in enumerate(trange[1:end-1])]])
NG = NoiseGrid(trange,Ws)

sol1, solend1 = @btime MSDE.sample(k_oop, u0, EM(false), NG)
sol2, solend2 = @btime MSDE.sample(k_oop, u0, MSDE.EulerMaruyama!(), Ws)
sol3, solend3 = @btime MSDE.sample(k_iip, u0, EM(false), NG)
sol4, solend4 = @btime MSDE.sample(k_iip2, u0, EM(false), NG)
sol5, solend5 = @btime MSDE.sample(k_iip2, u0, MSDE.EulerMaruyama!(), Ws)
sol6, solend6 = @btime MSDE.sample(k_iip2, u0, MSDE.EulerMaruyama!(), Ws, P=customP(p))

Z = collect(zip(1:length(trange), trange, Ws))
u = (1, 0.0, copy(u0))
dz = (0, 0.0, [0.0])
du = (0, 0.0, [0.0])
_, solend7 = MSDE.solve!(MSDE.EulerMaruyama!(), nothing, u, Z, customP(p))

@test solend1 ≈ solend2[3] rtol=1e-8
@test solend1 ≈ solend3 rtol=1e-8
@test solend1 ≈ solend4 rtol=1e-8
@test solend1 ≈ solend5[3] rtol=1e-8
@test solend1 ≈ solend6[3] rtol=1e-8
@test solend1 ≈ solend7[3] rtol=1e-8

u = (1, 0.0, copy(u0))
dz = (0, 0.0, [0.0])
@btime MSDE.solve!(MSDE.EulerMaruyama!(), nothing, deepcopy(u), Z, customP(p))
;


text:
using LinearAlgebra

function mc_gsa_estimate_total_indices(
    s::StochGalerkinODE,
    var_index,
    num_samples,
    interval_t,
    u0,
    tspan,
    p;
    kwargs...,
)
    @unpack ode_func!, dim_, vars = s
    A = zeros(num_samples, length(vars))
    B = zeros(num_samples, length(vars))
    for j in Base.OneTo(length(vars))
        for i in Base.OneTo(num_samples)
            A[i, j] = rand(vars[j])
            B[i, j] = rand(vars[j])
        end
    end
    C = deepcopy(A)
    C[:, var_index] .= B[:, var_index]
    Y_A = zeros(num_samples, length(interval_t), dim_)
    Y_C = zeros(num_samples, length(interval_t), dim_)
    for i in Base.OneTo(num_samples)
        sol_A = solve(
            ODEProblem(
                (du, u, p, t) -> ode_func!(du, u, p, (@view A[i, :]), t),
                u0,
                tspan,
                p,
            ),
            kwargs...,
        )
        sol_C = solve(
            ODEProblem(
                (du, u, p, t) -> ode_func!(du, u, p, (@view C[i, :]), t),
                u0,
                tspan,
                p,
            ),
            kwargs...,
        )
        for j in Base.OneTo(length(interval_t))
            Y_A[i, j, :] .= sol_A(interval_t[j])
            Y_C[i, j, :] .= sol_C(interval_t[j])
        end
    end
    sobol_indices = zeros(length(interval_t), dim_)
    for t_index in Base.OneTo(length(interval_t))
        for dim_index in Base.OneTo(dim_)
            Y_A_ = @view Y_A[:, t_index, dim_index]
            Y_C_ = @view Y_C[:, t_index, dim_index]
            expectation_not_xi = dot(Y_A_ .- Y_C_, Y_A_ .- Y_C_) / (2 * num_samples)
            variance = dot(Y_A_, Y_A_) / num_samples - (sum(Y_A_) / num_samples)^2
            sobol_indices[t_index, dim_index] = expectation_not_xi / variance
        end
    end
    return sobol_indices[2:end, :]
end

max_stars_repo_name: d-netto/PolyChaosODE

max_issues_repo_name: d-netto/PolyChaosODE

max_forks_repo_name: d-netto/PolyChaosODE


text:
using RandomCorrelationMatrices
using Base.Test

srand(12345)

println("Random correlation matrix")
n = 4
for η in [2,8,32,128]
    r = randcormatrix(n,η)
    @show η
    for i in 1:n
        for j in 1:n
            @printf("%5.2f ", r[i,j])
        end
        println()
    end

    ranges = Float64[]
    for s in 1:1000
        r = randcormatrix(n,η)
        r0 = r - diagm(diag(r))
        push!(ranges, maximum(vec(r0)) - minimum(vec(r0)))
    end
    @show mean(ranges)
    @show median(ranges)
end

println("\nRandom covariance matrix")
σ = [100, 200, 300, 400, 500]
@show σ
Σ = randcovmatrix(5, 10, σ)
L = chol(Σ,Val{:L})
samples = 100_000
results = zeros(5,samples)
for i in 1:samples
    results[:,i] = L*randn(5)
end
@show round(mean(results,2),1)
@show round(std(results,2),1)


text:
"""
High-order singular value decomposition (HO-SVD).
"""
function hosvd(tensor::AbstractArray{T,N}, core_dims::NTuple{N, Int};
               pad_zeros::Bool=false, compute_error::Bool=false) where {T,N}
    pad_zeros || _check_tensor(tensor, core_dims)

    factors = map(1:N) do i
        X = _col_unfold(tensor, i)
        f = eigen(Symmetric(X'X), max(1, size(X,2)-core_dims[i]+1):size(X,2)).vectors
        if pad_zeros && size(f, 2) < core_dims[i] # fill missing factors with zeros
            f = hcat(f, zeros(T, size(tensor, i), core_dims[i]-size(f, 2)))
        end
        mapslices(_check_sign, f, dims=1)
    end

    res = Tucker(factors, tensorcontractmatrices(tensor, factors))
    compute_error && _set_rel_residue(res, tensor)
    return res
end

hosvd(tensor::AbstractArray, r::Int;
      pad_zeros::Bool=false, compute_error::Bool=false) =
    hosvd(tensor, ntuple(_ -> r, ndims(tensor)),
          pad_zeros=pad_zeros, compute_error=compute_error)


text:
using Metatheory
using Metatheory.EGraphs

to_sexpr_pattern(p::PatLiteral) = "$(p.val)"
to_sexpr_pattern(p::PatVar) = "?$(p.name)"
function to_sexpr_pattern(p::PatTerm)
  e1 = join([p.head; to_sexpr_pattern.(p.args)], ' ')
  "($e1)"
end

to_sexpr(e::Symbol) = e
to_sexpr(e::Int64) = e
to_sexpr(e::Expr) = "($(join(to_sexpr.(e.args),' ')))"

function eggify(rules)
  egg_rules = []
  for rule in rules
    l = to_sexpr_pattern(rule.left)
    r = to_sexpr_pattern(rule.right)
    if rule isa SymbolicRule
      push!(egg_rules, "\tvec![rw!( \"$(rule.left) => $(rule.right)\" ; \"$l\" => \"$r\" )]")
    elseif rule isa EqualityRule
      push!(egg_rules, "\trw!( \"$(rule.left) == $(rule.right)\" ; \"$l\" <=> \"$r\" )")
    else
      println("Unsupported Rewrite Mode")
      @assert false
    end

  end
  return join(egg_rules, ",\n")
end

function rust_code(theory, query, params = SaturationParams())
  """
  use egg::{*, rewrite as rw};
  //use std::time::Duration;
  fn main() {
      let rules : &[Rewrite<SymbolLang, ()>] = &vec![
      $(eggify(theory))
      ].concat();

      let start = "$(to_sexpr(cleanast(query)))".parse().unwrap();
      let runner = Runner::default().with_expr(&start)
          // More options here https://docs.rs/egg/0.6.0/egg/struct.Runner.html
          .with_iter_limit($(params.timeout))
          .with_node_limit($(params.enodelimit))
          .run(rules);
      runner.print_report();
      let mut extractor = Extractor::new(&runner.egraph, AstSize);
      let (best_cost, best_expr) = extractor.find_best(runner.roots[0]);
      println!("best cost: {}, best expr {}", best_cost, best_expr);
  }
  """
end


text:
"""
    show

    Specialize Base.show(io::IO, question::Question) on this by passing a fixed header
"""
function show(io::IO, question::Question, header::Array{String,1})
    condition = ifelse(is_numeric(question.value), ">=", "==")
    colname = header[question.column]
    print(io, "Is $colname $condition $(question.value)?")
end

"""World's most elegant tree printing function."""
function print_tree(node::Leaf, spacing="")
    # Base case: we've reached a leaf
    println(spacing * "Predict" * repr(node.predictions))
end

function print_tree(node::Decision, spacing="")
    # Print the question at this node
    println(spacing * repr(node.question))

    # Call this function recursively on the true branch
    println(spacing * "--> True:")
    print_tree(node.true_branch, spacing * "  ")

    # Call this function recursively on the false branch
    println(spacing * "--> False:")
    print_tree(node.false_branch, spacing * "  ")
end