text:
addprocs(5)
################
#My local machine has 6 processors
#CYCLOPS is written to use them
#################
@everywhere basedir=homedir();
@everywhere netdir=string(basedir,"/Desktop/CYCLOPS_UPDATE_SHARE");
@everywhere cd(netdir);

using StatsBase
using MultivariateStats
using Distributions

@everywhere include("CYCLOPS_v6_2a_AutoEncoderModule_multi.jl")
@everywhere include("CYCLOPS_v6_2a_seed.jl")
@everywhere include("CYCLOPS_v6_2a_PreNPostprocessModule.jl")
@everywhere include("CYCLOPS_v6_2a_CircularStats_U.jl")
@everywhere include("CYCLOPS_v6_2a_MultiCoreModule_Smooth.jl")

using CYCLOPS_v6_2a_AutoEncoderModule_multi
using CYCLOPS_v6_2a_Seed
using CYCLOPS_v6_2a_PreNPostprocessModule
using CYCLOPS_v6_2a_CircularStats_U
using CYCLOPS_v6_2a_MultiCoreModule_Smooth

using PyPlot
#indirlaval=string(basedir,"/Documents/LungWork_March2015/LAVAL/UsableData");
#indirgrng=string(basedir,"/Documents/LungWork_March2015/GRNG/UsableData");
#homologuedir=string(basedir,"/Google Drive/BHTC_Homologues")

indirlaval=netdir;
indirgrng=netdir;
homologuedir=netdir;

############################################################
cd(indirlaval);
fullnonseed_data_laval=readcsv("AnnotatedLAVALData_March3_2015.csv");

cd(indirgrng);
fullnonseed_data_grng=readcsv("AnnotatedGRNGData_March3_2015.csv");

cd(homologuedir)
seed_homologues=readcsv("LungCyclerHomologues.csv");
homologue_symbol_list=seed_homologues[2:end,2];

############################################################
Frac_Var=0.85 # Set Number of Dimensions of SVD to maintain this fraction of variance
DFrac_Var=0.03 # Set Number of Dimensions of SVD to so that incremetal fraction of variance of var is at least this much
N_trials =20  # Number of random initial conditions to try for each optimization
MaxSeeds = 10000

total_background_num=20; # Number of background runs for global background refrence statistics
n_cores=5; # Number of machine cores
############################################################

alldata_data_l=fullnonseed_data_laval[2:end , 4:end];
alldata_data_l=Array{Float64}(alldata_data_l);

alldata_data_g=fullnonseed_data_grng[2:end , 4:end];
alldata_data_g=Array{Float64}(alldata_data_g);

alldata_symbols_g=fullnonseed_data_grng[2:end , 2];
nprobes=size(alldata_data_l)[1]
cutrank=nprobes-MaxSeeds

Seed_MinMean_l= (sort(vec(mean(alldata_data_l,2))))[cutrank]
Seed_MinMean_g= (sort(vec(mean(alldata_data_g,2))))[cutrank]

srand(123456);

seed_symbols_l, seed_data_l      =   getseed(fullnonseed_data_laval,homologue_symbol_list,.7,.14,Seed_MinMean_l,.975);
seed_data_dispersion_l           =   dispersion!(seed_data_l)
outs_l, norm_seed_data_l			   =	  GetEigenGenes(seed_data_dispersion_l,Frac_Var,DFrac_Var,30)

seed_symbols_g, seed_data_g      =   getseed(fullnonseed_data_grng,homologue_symbol_list,.7,.14,Seed_MinMean_g,.975);
seed_data_dispersion_g           =   dispersion!(seed_data_g)
outs_g, norm_seed_data_g				 =	 GetEigenGenes(seed_data_dispersion_g,Frac_Var,DFrac_Var,30)

#################
a_g =CYCLOPS_Order_multicore(outs_g,norm_seed_data_g,N_trials);
a_l =CYCLOPS_Order_multicore(outs_l,norm_seed_data_l,N_trials);

#a_g =CYCLOPS_Order(outs_g,norm_seed_data_g,N_trials);
#a_l =CYCLOPS_Order(outs_l,norm_seed_data_l,N_trials);

estimated_phaselist_g,bestnet_g,global_var_metrics_g=a_g;
estimated_phaselist_l,bestnet_l,global_var_metrics_l=a_l;

global_smooth_metrics_g            = smoothness_measures(seed_data_dispersion_g,norm_seed_data_g,estimated_phaselist_g)
global_metrics_g                   = global_var_metrics_g

global_smooth_metrics_l            = smoothness_measures(seed_data_dispersion_l,norm_seed_data_l,estimated_phaselist_l)
global_metrics_l                   = global_var_metrics_l

#pvals_g=multicore_backgroundstatistics_global_eigen(seed_data_dispersion_g,outs_g,N_trials,total_background_num,global_metrics_g)
#pvals_l=multicore_backgroundstatistics_global_eigen(seed_data_dispersion_l,outs_l,N_trials,total_background_num,global_metrics_l)

outdir=string(basedir,"/Google Drive/CYCLOPS_OUTPUT_FINAL_DISP/Lung/Nov2016_r1")
cd(outdir);

############################################################
function Filter_Cosinor_Output(cosdata::Array{Any,2},pval,rsq,ptt)
    significant_data=cosdata[append!([1],1 + findin(cosdata[2:end,5] .< pval,true)),:];
    phys_sig_data=significant_data[append!([1],1 + findin(significant_data[2:end,11] .> ptt,true)),:];
    strong_data=phys_sig_data[append!([1],1 + findin(phys_sig_data[2:end,10].> rsq,true)),:];
    strong_data
end
############################################################

estimated_phaselist_g=mod.(estimated_phaselist_g,2*pi)
estimated_phaselist_l=mod.(estimated_phaselist_l,2*pi)

cosinor_grng=Compile_MultiCore_Cosinor_Statistics(fullnonseed_data_grng,estimated_phaselist_g,4,24)
cosinor_laval=Compile_MultiCore_Cosinor_Statistics(fullnonseed_data_laval,estimated_phaselist_l,4,24)

sig_grng=Filter_Cosinor_Output(cosinor_grng,.05,0,1.66)
sig_laval=Filter_Cosinor_Output(cosinor_laval,.05,0,1.66)

common_g_l=intersect(sig_grng[2:end,1],sig_laval[2:end,1])
rows_laval_g_l=findin(sig_laval[:,1],common_g_l)
rows_grng_g_l=findin(sig_grng[:,1],common_g_l)

use_grng_g_l=mod.(sig_grng[rows_grng_g_l,6],2*pi)
use_laval_g_l=mod.(sig_laval[rows_laval_g_l,6],2*pi)

use_laval_adj1,estimated_phaselist_l_adj1=best_shift_cos2(use_laval_g_l,use_grng_g_l,estimated_phaselist_l,"radians")

scatter(use_laval_adj1,use_grng_g_l,s=.5,color="DarkBlue")

cosinor_laval=Compile_MultiCore_Cosinor_Statistics(fullnonseed_data_laval,estimated_phaselist_l_adj1,4,24)
####################################
# Align with Ebox (Par-BZip) phase
#####################################
eboxgenes=["DBP","HLF","TEF"]

eboxphases_laval=cosinor_laval[findin(cosinor_laval[:,2],eboxgenes),:]
eboxphases_grng=cosinor_grng[findin(cosinor_grng[:,2],eboxgenes),:]

laval_criteria=((eboxphases_laval[:,4].<.05) & ( eboxphases_laval[:,11].>1.25) & ( eboxphases_laval[:,9].>100))
grng_criteria=((eboxphases_grng[:,4].<.05) & ( eboxphases_grng[:,11].>1.25) & ( eboxphases_laval[:,9].>100))

eboxphases_laval=eboxphases_laval[findin(laval_criteria,true),:]
eboxphases_grng=eboxphases_grng[findin(grng_criteria,true),:]

eboxphases_laval=Array{Float64}(eboxphases_laval[:,6])
eboxphases_grng=Array{Float64}(eboxphases_grng[:,6])

eboxphases=append!(eboxphases_laval,eboxphases_grng)
eboxphase=Circular_Mean(eboxphases)

estimated_phaselist_g_adj_final=mod.(estimated_phaselist_g .- eboxphase+(pi),2*pi)  ##ebox genes in mouse lung peak on average CT 11.5
estimated_phaselist_l_adj_final=mod.(estimated_phaselist_l_adj1 .- eboxphase+(pi),2*pi)  ##ebox genes in mouse lung peak on average CT 11.5

Nday=length(findin((estimated_phaselist_l_adj_final .>pi),true))
Nnight=length(findin((estimated_phaselist_l_adj_final .<pi),true))

if (Nday<Nnight)
    estimated_phaselist_l_adj_final=mod(2*pi-estimated_phaselist_l_adj_final,2*pi)
    estimated_phaselist_g_adj_final=mod(2*pi-estimated_phaselist_g_adj_final,2*pi)
end

###################################
cosinor_grng=Compile_MultiCore_Cosinor_Statistics(fullnonseed_data_grng,estimated_phaselist_g_adj_final,4,24)
cosinor_laval=Compile_MultiCore_Cosinor_Statistics(fullnonseed_data_laval,estimated_phaselist_l_adj_final,4,24)

nsig_grng=Filter_Cosinor_Output(cosinor_grng,.05,0,1.66)
nsig_laval=Filter_Cosinor_Output(cosinor_laval,.05,0,1.66)

common_g_l=intersect(nsig_grng[2:end,1],nsig_laval[2:end,1])

cd(outdir);
#writecsv("Laval_Sample_Phaselist.csv",estimated_phaselist_l_adj_final);
#writecsv("GRNG_Sample_Phaselist.csv",estimated_phaselist_g_adj_final);
#writecsv("Laval_Cosinor_Output.csv",cosinor_laval);
#writecsv("GRNG_Cosinor_Output.csv",cosinor_grng);

############################################################
sig_laval=Filter_Cosinor_Output(cosinor_laval,.05,0,1.666)
sig_grng=Filter_Cosinor_Output(cosinor_grng,.05,0,1.666)

############################################################
common_g_l=intersect(sig_grng[2:end,1],sig_laval[2:end,1])

rows_laval_g_l=findin(sig_laval[:,1],common_g_l)
rows_grng_g_l=findin(sig_grng[:,1],common_g_l)

use_grng_g_l=sig_grng[rows_grng_g_l,6]
use_laval_g_l=sig_laval[rows_laval_g_l,6]
############################################################
# Acrophase Comparison
############################################################
close()
ylabp=[0,pi/2,pi,3*pi/2,2*pi]
ylabs=[0, "","π","","2π"]
xlabp=[0,pi/2,pi,3*pi/2,2*pi]
xlabs=[0, "","π","","2π"]

scatter(mod(use_grng_g_l,2*pi),mod(use_laval_g_l,2*pi),s=.5,color="DarkBlue")
xticks(xlabp, xlabs, fontsize=18)
yticks(ylabp, ylabs, fontsize=18)
#xlabel("GRNG Acrophase (radians)", fontsize=18)
#ylabel("Laval Acrophase (radians)", fontsize=18)

#scatter(use_grng_g_l,use_laval_g_l,s=.5,color="DarkBlue")

############################################################
# Sample Historgrams
############################################################
close()
subplot(1,2,1)
plt[:hist](estimated_phaselist_g_adj_final,normed=1,alpha=.5,color="DarkGreen")
xlabel("Phase of Sample Collection (GRNG)", fontsize=18)
ylabel("Probability", fontsize=18)
subplot(1,2,2)
plt[:hist](estimated_phaselist_l_adj_final,normed=1,alpha=.5,color="DarkBlue")
ylabel("Probability", fontsize=18)
xlabel("Phase of Sample Collection (Laval)", fontsize=18)
############################################################
close()
plt[:hist](estimated_phaselist_l_adj_final,normed=1,alpha=.5,color="DarkBlue")
xticks(xlabp, xlabs)
ylabel("Probability", fontsize=18)
xlabel("Phase of Sample Collection (Laval)", fontsize=18)
############################################################
# Sample Historgrams - Laval - No annotations
############################################################
close()
plt[:hist](estimated_phaselist_l_adj_final,normed=1,alpha=.5,color="DarkBlue")
xticks(xlabp, xlabs)
############################################################
# Sample Historgrams - All- No annotations
############################################################
close()
plt[:hist](vcat(estimated_phaselist_l_adj_final,estimated_phaselist_g_adj_final),normed=1,alpha=.5,color="DarkBlue")
xticks(xlabp, xlabs)
ylabel("Probability", fontsize=18)
xlabel("Phase of Sample Collection (Laval & GRNG)", fontsize=18)
###################################
r_grng=Filter_Cosinor_Output(cosinor_grng,.05,.65,1.66)
r_laval=Filter_Cosinor_Output(cosinor_laval,.05,.65,1.66)

common_g_l_r=intersect(r_grng[2:end,2],r_laval[2:end,2])
common_g_l_r[common_g_l_r .!=""]

common_g_l_r=intersect(r_grng[2:end,1],r_laval[2:end,1])

clockgenes=["CLOCK","CRY2","CRY1","DBP","TEF","HLF","EFNB2","ADAM9","E4F1","ACE","NR1D1","NR1D2","RORA"]
aclockrows=findin(alldata_symbols_g,clockgenes);
clockrows=aclockrows
clockrows=aclockrows[[23,4,27,21,6,16,8,20,17]]
clockrows=aclockrows[[23,4,27,8,20,17]]
clockdata_g=alldata_data_g[clockrows,:];
clockdata_l=alldata_data_l[clockrows,:];

cosinor_clock_g=cosinor_grng[clockrows+1,:];  ## the +1 is because the first row of this array is the headings
cosinor_clock_l=cosinor_laval[clockrows+1,:];  ## the +1 is because the first row of this array is the headings

clockannotations=alldata_symbols_g[clockrows,:];

clockdata_g=alldata_data_g[clockrows,:];
clockdata_l=alldata_data_l[clockrows,:];

cosinor_clock_g=cosinor_grng[clockrows+1,:];  ## the +1 is because the first row of this array is the headings
cosinor_clock_l=cosinor_laval[clockrows+1,:];  ## the +1 is because the first row of this array is the headings

clockannotations=alldata_symbols_g[clockrows,:];
close()

close()

xlabp=[0,pi/2,pi,3*pi/2,2*pi]
xlabs=[0, "","π","","2π"]

for n in 1:6
          m1=2*n-1
          m2=2*n
          subplot(2,6,m1)

          ymean=mean(clockdata_g[n,:])
          ymax=minimum([ymean+3*std(clockdata_g[n,:]),maximum(clockdata_g[n,:])])
          ymin=maximum([0,ymean-3*std(clockdata_g[n,:])])

          ymaxstar=Integer(100*div(ymax,100))
          yminstar=Integer(100*div(ymin,100))
          ymedstar=Integer((yminstar+ymaxstar)/2)

          ylabp=[yminstar,ymedstar,ymaxstar]
          ylabs=[yminstar,ymedstar,ymaxstar]

          axis([0,2*pi,ymin,ymax])
          xticks(xlabp, xlabs)
          yticks(ylabp, ylabs, fontsize=10)

          scatter(estimated_phaselist_g_adj_final,clockdata_g[n,:],alpha=.75,s=7,color="DarkBlue")

          title(clockannotations[n], fontsize=18)
          xticks(xlabp, xlabs)
          yticks(ylabp, ylabs, fontsize=10)

          PrbPhase,PrbAmp,PrbMean=cosinor_clock_g[n,[6,7,8]]
          sest=linspace(0,2*pi,100)
          synth=PrbAmp*cos.(sest-PrbPhase)+PrbMean
          plot(sest,synth,"r-",lw=2)

          subplot(2,6,m2)

          ymean=0
          ymax=0
          ymin=0

          ymean=mean(clockdata_l[n,:])
          ymax=minimum([ymean+3*std(clockdata_l[n,:]),maximum(clockdata_l[n,:])])
          ymin=maximum([0,ymean-3*std(clockdata_l[n,:])])
          println(yminstar,"   ",ymedstar,"    ",ymaxstar)

          ymaxstar=Integer(100*div(ymax,100))
          yminstar=Integer(100*div(ymin,100))
          ymedstar=Integer((yminstar+ymaxstar)/2)

          ylabp=[yminstar,ymedstar,ymaxstar]
          ylabs=[yminstar,ymedstar,ymaxstar]

          axis([0,2*pi,ymin,ymax])
          xticks(xlabp, xlabs)
          yticks(ylabp, ylabs, fontsize=10)

          scatter(estimated_phaselist_l_adj_final,clockdata_l[n,:],alpha=.75,s=7,color="DarkGreen")

          title(clockannotations[n], fontsize=18)

          PrbPhase,PrbAmp,PrbMean=cosinor_clock_l[n,[6,7,8]]
          sest=linspace(0,2*pi,100)
          synth=PrbAmp*cos.(sest-PrbPhase)+PrbMean
          plot(sest,synth,"r-",lw=2)

end

#suptitle("",fontsize=22)

clockgenes=["KDR","ACE","MAP4K3","SMAD7","AQP4"]
clockrows=findin(alldata_symbols_g,clockgenes);
clockrows=clockrows[[1,4,5,7,11]]
clockdata_g=alldata_data_g[clockrows,:];
clockdata_l=alldata_data_l[clockrows,:];

cosinor_clock_g=cosinor_grng[clockrows+1,:];  ## the +1 is because the first row of this array is the headings
cosinor_clock_l=cosinor_laval[clockrows+1,:];  ## the +1 is because the first row of this array is the headings

clockannotations=alldata_symbols_g[clockrows,:];
close()
#############################

for n in 1:5
          m1=2*n-1
          m2=2*n
          subplot(5,2,m1)

          ymean=mean(clockdata_g[n,:])
          ymax=minimum([ymean+3*std(clockdata_g[n,:]),maximum(clockdata_g[n,:])])
          ymin=maximum([0,ymean-3*std(clockdata_g[n,:])])

          ymaxstar=Integer(100*div(ymax,100))
          yminstar=Integer(100*div(ymin,100))
          ymedstar=Integer((yminstar+ymaxstar)/2)

          ylabp=[yminstar,ymedstar,ymaxstar]
          ylabs=[yminstar,ymedstar,ymaxstar]

          axis([0,2*pi,ymin,ymax])
          xticks(xlabp, xlabs)
          yticks(ylabp, ylabs, fontsize=10)

          scatter(estimated_phaselist_g_adj_final,clockdata_g[n,:],alpha=.75,s=7,color="DarkBlue")

          title(clockannotations[n], fontsize=18)
          xticks(xlabp, xlabs)
          yticks(ylabp, ylabs, fontsize=10)

          PrbPhase,PrbAmp,PrbMean=cosinor_clock_g[n,[6,7,8]]
          sest=linspace(0,2*pi,100)
          synth=PrbAmp*cos.(sest-PrbPhase)+PrbMean
          plot(sest,synth,"r-",lw=2)

          subplot(5,2,m2)

          ymean=0
          ymax=0
          ymin=0

          ymean=mean(clockdata_l[n,:])
          ymax=minimum([ymean+3*std(clockdata_l[n,:]),maximum(clockdata_l[n,:])])
          ymin=maximum([0,ymean-3*std(clockdata_l[n,:])])
          println(yminstar,"   ",ymedstar,"    ",ymaxstar)

          ymaxstar=Integer(100*div(ymax,100))
          yminstar=Integer(100*div(ymin,100))
          ymedstar=Integer((yminstar+ymaxstar)/2)

          ylabp=[yminstar,ymedstar,ymaxstar]
          ylabs=[yminstar,ymedstar,ymaxstar]

          axis([0,2*pi,ymin,ymax])
          xticks(xlabp, xlabs)
          yticks(ylabp, ylabs, fontsize=10)

          scatter(estimated_phaselist_l_adj_final,clockdata_l[n,:],alpha=.75,s=7,color="DarkGreen")

          title(clockannotations[n], fontsize=18)

          PrbPhase,PrbAmp,PrbMean=cosinor_clock_l[n,[6,7,8]]
          sest=linspace(0,2*pi,100)
          synth=PrbAmp*cos.(sest-PrbPhase)+PrbMean
          plot(sest,synth,"r-",lw=2)

end
########################################
###########################
# Demonstration that allowing more itterations of neural network does not much change fit
#########################
function refine_order(l_outs::Integer,l_norm_seed_data::Array{Float64,2},NET,epochs=50000)
    scalefactor=1;
    besterror=10E20
    bestnet=NET;

    Train_Momentum_Stochastic!(l_norm_seed_data,NET,10,.1,.5,0.0005,epochs)
    global_circ_sse=Train_Bold!(l_norm_seed_data,NET,.1,0.0005,epochs)
    besterror=global_circ_sse
    print(global_circ_sse)

    total_sse,onedlinear_sse=Default_Metrics(l_norm_seed_data,l_outs)
    besterror=besterror*2
    global_metrics=[besterror, besterror/(total_sse-besterror), besterror/(onedlinear_sse-besterror)]

    estimated_phaselist=ExtractPhase(l_norm_seed_data,bestnet);
    estimated_phaselist=mod.(estimated_phaselist + 2*pi,2*pi)
    estimated_phaselist,bestnet,global_metrics
end
########################
estimated_phaselist_gh,bestnet_gh,global_var_metrics_gh=refine_order(outs_g,norm_seed_data_g,bestnet_g,30000);
estimated_phaselist_lh,bestnet_lh,global_var_metrics_lh=refine_order(outs_l,norm_seed_data_l,bestnet_l,30000);

close()
subplot(1,2,1)
scatter(estimated_phaselist_g,estimated_phaselist_gh)
subplot(1,2,2)
scatter(estimated_phaselist_l,estimated_phaselist_lh)

end

max_stars_repo_name: AdithyaVenkateshMohan/CYCLOPS

max_issues_repo_name: AdithyaVenkateshMohan/CYCLOPS

max_forks_repo_name: AdithyaVenkateshMohan/CYCLOPS


text:
"""
This file contains the code for the work by

P. S. Thomas, G. Theocharous, and M. Ghavamzadeh.
High Confidence Off-Policy Evaluation.
In Proceedings of the Twenty-Ninth AAAI Conference on Artificial Intelligence, 2015.

Code repo: https://github.com/yashchandak/UnO
"""

using Statistics
using Random
using Bootstrap

function OPE_PDIS(data, props, ξr=0.0)
    """
    Off policy evaluation using per-decision importance sampling

    ----
    Args

    data: Trajectories x Steps x (Importance ratio, Reward)
    prop: Meta properties of the data
    ξr: can be used to shift the rewards such that it is always non-positive or non-negative
    """

    # Extract the meta properties of the data
    n, horizon, d = size(data)
    γ = props["gamma"]

    total = fill(0.0, n)                        # Stores off-policy estimates for each trajectory
    for i in 1:1:n
        offreturn = 0.0                         # Initialize off-policy estimate for a trajectory
        ρ = 1.0                                 # Inititalize the importance ratio
        for j in 1:1:horizon
            ρ = ρ*data[i,j,1]                   # Product of past IS ratios
            reward = data[i,j,2] - ξr
            offreturn += γ^(j-1) * ρ * reward   # j-1 because first power should be 0
        end;
        total[i] = offreturn
    end;
    return total
end;

function OPE_IS(data, props, ξG=0.0)
    """
    Off policy evaluation using full-trajectory importance sampling

    ----
    Args

    data: Trajectories x Steps x (Importance ratio, Reward)
    prop: Meta properties of the data
    ξG: can be used to shift the returns such that it is always non-positive or non-negative
    """

    # Extract the meta properties of the data
    n, horizon, d = size(data)
    γ = props["gamma"]

    total = fill(0.0, n)
    for i in 1:1:n
        ρ = 1.0
        returns = - ξG
        for j in 1:1:horizon
            # Full trajectory IS
            ρ = ρ*data[i,j,1]
            returns += γ^(j-1) * data[i,j,2]
        end;
        total[i] = ρ * returns
    end;
    return total
end;

function CI_lower(data, props, δ=0.05)
    """
    Lower bound on the mean using concentration inequalities

    Args
    ----

    data: Trajectories x Steps x (Importance ratio, Reward)
    prop: Meta properties of the data
    δ: Acceptable failure rate
    """

    # Extract meta-properties
    n, horizon, d = size(data)
    γ, T = props["gamma"], props["horizon"]
    rmin, gmin = props["rmin"], props["Gmin"]
    geom_sum = (1 - γ^T)/(1 - γ)

    # Get the importance sampled returns ** with control variate = rmin **
    estimates = OPE_PDIS(data, props, rmin)
    if any(x -> x<0, estimates)
        println("Bound will not work, items in array less than 0 (CI Lower)")
    end;

    #########################################################
    # Lower bound proposed by Thomas et al.  (2015)

    function lower_bound(D, c, n_post)
        Y = clamp.(D, 0.0, c)
        term1 = mean(Y)
        term2 = (7*c*log(2/δ)) / (3*(n_post - 1))
        term3 = sqrt( (2 * var(Y) * log(2/δ) )/ n_post) # Sample variance form from Maurer(2009)

        lb = term1 - term2 - term3                      # Thm1 from Thomas et al. (2015)
        lb = lb + rmin * geom_sum                       # add back control to the bound (Thm 5 from Chandak et al. (2021))
        return lb
    end;

    #########################################################
    # Prepare data splits into dpre and dpost
    # Size based on the split-value suggested by Thomas et al.  (2015)

    n_pre = Int(floor(0.05 * n))
    n_post = n - n_pre

    perm_idxs = shuffle(collect(1:n))
    D_pre = estimates[perm_idxs[1:n_pre]]
    D_post = estimates[perm_idxs[n_pre + 1:n]]

    ##########################################################
    # Use D_pre to search for best truncation point
    # using basic random search in 1D

    min_obs, max_obs = minimum(D_pre), maximum(D_pre)
    cs = (rand(100) .* (max_obs - min_obs)) .+ min_obs
    v_max, c_max = -Inf, cs[1]

    # Search for c that maximizes the lower bound
    for c in cs
        # Notice the use of n_post instead of n_pre for the size of data.
        # See the work by Thomas et al.  (2015) for discussion.
        val = lower_bound(D_pre, c, n_post)
        if val > v_max
            v_max, c_max = val, c
        end;
    end;

    ########################################################
    # Use D_post to obtain the final bound by Thomas et al. (2015)
    lb = lower_bound(D_post, c_max, n_post)

    # CI lower cannot be less than minimum possible returns
    lb = max(lb, geom_sum*rmin, gmin)

    return lb
end;

function CI_upper(data, props, δ=0.05)
    """
    Upper bound on the mean using concentration inequalities

    Args
    ----

    data: Trajectories x Steps x (Importance ratio, Reward)
    prop: Meta properties of the data
    δ: Acceptable failure rate
    """

    # Extract meta-properties
    n, horizon, d = size(data)
    γ, T = props["gamma"], props["horizon"]
    rmax, gmax = props["rmax"], props["Gmax"]
    geom_sum = (1 - γ^T)/(1 - γ)

    # Get the importance sampled returns ** with control variate **
    estimates = OPE_PDIS(data, props, rmax)
    if any(x -> x>0, estimates)
        println("Bound will not work, items in array more than 0 (CI Upper)")
    end;

    #########################################################
    # Upper bound extension for the bound by Thomas et al.  (2015)
    # See Thm 8 in the work by Chandak et al. (2021) for details.

    function upper_bound(D, c, n_post)
        Y = clamp.(D, c, 0.0)
        term1 = mean(Y)
        term2 = (7*abs(c)*log(2/δ)) / (3*(n_post - 1))
        term3 = sqrt( (2 * var(Y) * log(2/δ) )/ n_post) # Sample variance form from Maurer(2009)
        # println(term1, term2, term3)
        ub = term1 + term2 + term3                      # Thm1 from Thomas et al.  (2015)
        ub = ub + rmax*geom_sum                         # add back control to the bound (Thm 5 by Chandak et al. (2021))
        return ub
    end;

    #########################################################
    # Prepare data splits into dpre and dpost
    # Size based on the split-value suggested by Thomas et al.  (2015)

    n_pre = Int(floor(0.05 * n))
    n_post = n - n_pre

    perm_idxs = shuffle(collect(1:n))
    D_pre = estimates[perm_idxs[1:n_pre]]
    D_post = estimates[perm_idxs[n_pre + 1:n]]

    ##########################################################
    # Use D_pre to search for best truncation point
    # using basic random search in 1D

    min_obs, max_obs = minimum(D_pre), maximum(D_pre)
    cs = (rand(100) .* (max_obs - min_obs)) .+ min_obs
    v_min, c_min = +Inf, cs[1]

    # Search for c that minimizes the upper bound
    for c in cs
        # Notice the use of n_post instead of n_pre for the size of data
        # See the work by Thomas et al.  (2015) for discussion.
        val = upper_bound(D_pre, c, n_post)
        if val < v_min
            v_min, c_min = val, c
        end;
    end;

    ########################################################
    # Use D_post to obtain the final bound
    ub = upper_bound(D_post, c_min, n_post)

    # CI should not be higher than the maximum possible returns
    ub = min(ub, rmax*geom_sum, gmax)

    return ub
end;

function HCOPE_Bootstrap(data, props, δ=0.05)
    """
    Upper and Lower (approx) bound on the mean using Bootstrap

    Args
    ----

    data: Trajectories x Steps x (Importance ratio, Reward)
    prop: Meta properties of the data
    δ: Acceptable failure rate
    """

    full_size, _, _ = size(data)

    IS_data = OPE_PDIS(data, props)
    bs = bootstrap(mean, IS_data, BasicSampling(100))   # Bootstrap mean's distribution
    # bci = confint(bs, BCaConfInt(0.95))               # Bootstrap CI on the mean
    bci = confint(bs, PercentileConfInt(1 - δ))         # Bootstrap CI on the mean

    estimate, lb, ub = bci[1]                           # Obtain the estiamte, lower bound, and the upper bound

    γ, T = props["gamma"], props["horizon"]
    rmin, gmin = props["rmin"], props["Gmin"]
    rmax, gmax = props["rmax"], props["Gmax"]
    geom_sum = (1 - γ^T)/(1 - γ)

    # Intersect bootstrap intervals with maximum and minimum possible return values
    lb = max(lb, gmin, geom_sum*rmin)
    ub = min(ub, gmax, geom_sum*rmax)

    return [estimate, lb, ub]
end;

ext: jl

max_stars_repo_name: yashchandak/UnO

max_issues_repo_name: yashchandak/UnO

max_forks_repo_name: yashchandak/UnO


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -51.708688423784224 <= q <= 53.49505517188359)

                     add_NL_constraint(m, :(softplus(0.12651985253467002 + 0.8800990177661969*softplus(0.9367289958343452 + 0.010085610540752743*softplus(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.46191185446397043*softplus(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.4430556996324988*softplus(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.9244167282970754*softplus(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.15008923933940022*softplus(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + 0.9059509477084742*softplus(-0.9076856947519341 + 0.6427398296015978*softplus(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + 0.5974696839876872*softplus(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + 0.3127996929073351*softplus(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.8921846648994118*softplus(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.8861141904898706*softplus(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + 0.7662733112940292*softplus(0.014219790982957559 + -0.1830460372609095*softplus(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + 0.18937048380167898*softplus(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.48856019606524415*softplus(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.6243584215460145*softplus(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + -0.6378733350620909*softplus(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + -0.3066544631255512*softplus(0.37873800689777815 + 0.3786798527526116*softplus(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.15214183164812667*softplus(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + 0.4531653277321981*softplus(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + -0.2698943947473773*softplus(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + -0.8205030207972928*softplus(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + 0.673637701918695*softplus(0.41028851710489267 + 0.3789873528571355*softplus(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.13395473418456794*softplus(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.9211828416113521*softplus(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.005265087561311166*softplus(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.6679558602332873*softplus(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4])))) + softplus(-0.8950204772394628 + -0.19603919659654778*softplus(0.9367289958343452 + 0.010085610540752743*softplus(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.46191185446397043*softplus(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.4430556996324988*softplus(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.9244167282970754*softplus(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.15008923933940022*softplus(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + -0.902347490498435*softplus(-0.9076856947519341 + 0.6427398296015978*softplus(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + 0.5974696839876872*softplus(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + 0.3127996929073351*softplus(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.8921846648994118*softplus(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.8861141904898706*softplus(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + -0.30440629677970144*softplus(0.014219790982957559 + -0.1830460372609095*softplus(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + 0.18937048380167898*softplus(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.48856019606524415*softplus(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.6243584215460145*softplus(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + -0.6378733350620909*softplus(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + 0.6770808003403048*softplus(0.37873800689777815 + 0.3786798527526116*softplus(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.15214183164812667*softplus(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + 0.4531653277321981*softplus(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + -0.2698943947473773*softplus(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + -0.8205030207972928*softplus(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + 0.03232239144278681*softplus(0.41028851710489267 + 0.3789873528571355*softplus(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.13395473418456794*softplus(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.9211828416113521*softplus(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.005265087561311166*softplus(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.6679558602332873*softplus(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4])))) + softplus(-0.11076390433121786 + 0.7661103402200116*softplus(0.9367289958343452 + 0.010085610540752743*softplus(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.46191185446397043*softplus(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.4430556996324988*softplus(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.9244167282970754*softplus(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.15008923933940022*softplus(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + -0.9009199905889229*softplus(-0.9076856947519341 + 0.6427398296015978*softplus(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + 0.5974696839876872*softplus(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + 0.3127996929073351*softplus(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.8921846648994118*softplus(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.8861141904898706*softplus(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + 0.5631507797475948*softplus(0.014219790982957559 + -0.1830460372609095*softplus(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + 0.18937048380167898*softplus(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.48856019606524415*softplus(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.6243584215460145*softplus(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + -0.6378733350620909*softplus(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + -0.4163802251888926*softplus(0.37873800689777815 + 0.3786798527526116*softplus(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.15214183164812667*softplus(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + 0.4531653277321981*softplus(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + -0.2698943947473773*softplus(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + -0.8205030207972928*softplus(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + -0.17662226846197893*softplus(0.41028851710489267 + 0.3789873528571355*softplus(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.13395473418456794*softplus(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.9211828416113521*softplus(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.005265087561311166*softplus(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.6679558602332873*softplus(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4])))) + softplus(0.8242960442017702 + -0.7506191216499181*softplus(0.9367289958343452 + 0.010085610540752743*softplus(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.46191185446397043*softplus(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.4430556996324988*softplus(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.9244167282970754*softplus(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.15008923933940022*softplus(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + 0.31521888864028513*softplus(-0.9076856947519341 + 0.6427398296015978*softplus(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + 0.5974696839876872*softplus(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + 0.3127996929073351*softplus(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.8921846648994118*softplus(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.8861141904898706*softplus(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + 0.5308637261158768*softplus(0.014219790982957559 + -0.1830460372609095*softplus(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + 0.18937048380167898*softplus(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.48856019606524415*softplus(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.6243584215460145*softplus(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + -0.6378733350620909*softplus(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + -0.6654841763336146*softplus(0.37873800689777815 + 0.3786798527526116*softplus(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.15214183164812667*softplus(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + 0.4531653277321981*softplus(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + -0.2698943947473773*softplus(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + -0.8205030207972928*softplus(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + 0.2645412007156045*softplus(0.41028851710489267 + 0.3789873528571355*softplus(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.13395473418456794*softplus(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.9211828416113521*softplus(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.005265087561311166*softplus(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.6679558602332873*softplus(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4])))) + softplus(-0.6657453597374303 + 0.46323784626870257*softplus(0.9367289958343452 + 0.010085610540752743*softplus(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.46191185446397043*softplus(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.4430556996324988*softplus(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.9244167282970754*softplus(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.15008923933940022*softplus(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + 0.7128966868012561*softplus(-0.9076856947519341 + 0.6427398296015978*softplus(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + 0.5974696839876872*softplus(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + 0.3127996929073351*softplus(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.8921846648994118*softplus(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.8861141904898706*softplus(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + 0.43240821263487783*softplus(0.014219790982957559 + -0.1830460372609095*softplus(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + 0.18937048380167898*softplus(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.48856019606524415*softplus(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.6243584215460145*softplus(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + -0.6378733350620909*softplus(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + -0.6998749584965918*softplus(0.37873800689777815 + 0.3786798527526116*softplus(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.15214183164812667*softplus(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + 0.4531653277321981*softplus(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + -0.2698943947473773*softplus(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + -0.8205030207972928*softplus(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4]))) + 0.6177298271304656*softplus(0.41028851710489267 + 0.3789873528571355*softplus(-0.5507415348587976 + -0.9382157052845588*$(x[1]) + 0.32792621601464056*$(x[2]) + 0.04402367930308282*$(x[3]) + 0.8348925442898456*$(x[4])) + -0.13395473418456794*softplus(0.4569151270063867 + 0.9855428533575652*$(x[1]) + 0.27576336558204195*$(x[2]) + 0.024724286249658167*$(x[3]) + -0.4115946180479595*$(x[4])) + -0.9211828416113521*softplus(0.5124722980797523 + -0.021522352811447654*$(x[1]) + 0.09961216478582813*$(x[2]) + -0.5936461652265081*$(x[3]) + 0.7973439553626904*$(x[4])) + 0.005265087561311166*softplus(0.16894693163358276 + -0.35673489881937037*$(x[1]) + 0.5082818593140526*$(x[2]) + 0.7581208342644015*$(x[3]) + -0.6603822729791307*$(x[4])) + 0.6679558602332873*softplus(-0.5086292355642796 + -0.9220327700183186*$(x[1]) + 0.16243082406428888*$(x[2]) + 0.7781879917830965*$(x[3]) + 0.45589624224566805*$(x[4])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using Test, SummationByPartsOperators
using LinearAlgebra

function accuracy_test!(res, ufunc, dufunc, D)
    u = compute_coefficients(ufunc, D)
    du = compute_coefficients(dufunc, D)
    mul!(res, D, u)
    maximum(abs, du-res) < 5*length(res)*eps(eltype(res))
end

# Accuracy Tests
for T in (Float32, Float64)
    xmin = -one(T)
    xmax = one(T)

    for N in 2 .^ (3:6)
        D = fourier_derivative_operator(xmin, xmax, N)
        println(devnull, D)
        @test SummationByPartsOperators.derivative_order(D) == 1
        @test issymmetric(D) == false
        u = compute_coefficients(zero, D)
        res = D*u
        for k in 0:(N÷2)-1
            ufunc = x->sinpi(k*x)
            dufunc = x->typeof(x)(k*π)*cospi(k*x)
            @test accuracy_test!(res, ufunc, dufunc, D)
            xplot, duplot = evaluate_coefficients(res, D)
            @test maximum(abs, duplot - dufunc.(xplot)) < 5N*eps(T)
            @test abs(integrate(u, D)) < N*eps(T)

            ufunc = x->cospi(k*x)
            dufunc = x->-typeof(x)(k*π)*sinpi(k*x)
            @test accuracy_test!(res, ufunc, dufunc, D)
            xplot, duplot = evaluate_coefficients(res, D)
            @test maximum(abs, duplot - dufunc.(xplot)) < 5N*eps(T)
            @test abs(integrate(u, D)) < N*eps(T)
        end
    end
end

# (Super) Spectral Viscosity
source_SV = (Tadmor1989(), MadayTadmor1989(), TadmorWaagan2012Standard(), TadmorWaagan2012Convergent())
source_SSV = (Tadmor1993(),)

for T in (Float32, Float64), source in source_SV
    xmin = -one(T)
    xmax = one(T)

    for N in 2 .^ (3:6)
        D = fourier_derivative_operator(xmin, xmax, N)
        println(devnull, D)
        @test issymmetric(D) == false

        Di = dissipation_operator(source, D)
        println(devnull, Di)
        @test issymmetric(Di) == true
        Di_full = Matrix(Di)
        @test maximum(abs, Di_full-Di_full') < 80*eps(T)

        @test maximum(eigvals(Symmetric(Di_full))) < 10N*eps(T)
    end
end

for T in (Float32, Float64), source in source_SSV
    xmin = -one(T)
    xmax = one(T)

    for N in 2 .^ (3:6), order in 1:3
        D = fourier_derivative_operator(xmin, xmax, N)
        println(devnull, D)
        @test issymmetric(D) == false

        Di = dissipation_operator(source, D, order=order)
        println(devnull, Di)
        @test issymmetric(Di) == true
        Di_full = Matrix(Di)
        @test maximum(abs, Di_full-Di_full') < 80*eps(T)

        @test maximum(eigvals(Symmetric(Di_full))) < 15N*eps(T)
    end
end

# Modal Filtering
for T in (Float32, Float64), filter_type in (ExponentialFilter(),)
    xmin = -one(T)
    xmax = one(T)

    for N in 2 .^ (1:4)
        D = fourier_derivative_operator(xmin, xmax, N)
        filter! = ConstantFilter(D, filter_type)
        u = compute_coefficients(zero, D)
        res = D*u
        for k in 1:N-1
            compute_coefficients!(u, x->exp(sinpi(x)), D)
            norm2_u = integrate(u->u^2, u, D)
            filter!(u)
            @test integrate(u->u^2, u, D) <= norm2_u
        end
    end
end


text:
using DrWatson
quickactivate(@__DIR__)
@show projectdir()
##
using RCall
using CSV
using BSON
using DataFrames
using Dates
using Plots
using Underscores
using StatsBase
using ColorSchemes
plotlyjs()
##

contacts = BSON.load( projectdir("data/contacts/int/mean_contact_rates_aggregation=14.bson") )
 = contacts["Denmark"].date
diff(rw_step)|>unique

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
module AttachConstraintsTests

using Test
using FillArrays
using LinearAlgebra
using Gridap.Helpers
using Gridap.Fields
using Gridap.ReferenceFEs
using Gridap.Arrays
using Gridap.Integration
using Gridap.CellData
using Gridap.TensorValues

#using Gridap.Fields: MockField, MockBasis, OtherMockBasis

ncells = 10
ndofs = 3
ndofs_c = 4
cellmat = [ rand(ndofs,ndofs) for cell in 1:ncells]
cellvec = [ rand(ndofs) for cell in 1:ncells]
cellmatvec = pair_arrays(cellmat,cellvec)
cellconstr = [ rand(ndofs_c,ndofs) for cell in 1:ncells]

a = attach_constraints_rows(cellvec,cellconstr)
r = map( (vec,constr) -> constr*vec ,cellvec,cellconstr)
test_array(a,r)

a = attach_constraints_rows(cellmat,cellconstr)
r = map( (mat,constr) -> constr*mat ,cellmat,cellconstr)
test_array(a,r)

a = attach_constraints_rows(cellmatvec,cellconstr)
r = map( (mat,vec,constr) -> (constr*mat,constr*vec),cellmat,cellvec,cellconstr)
test_array(a,r)

a = attach_constraints_cols(cellmat,cellconstr)
r = map( (mat,constr) -> mat*transpose(constr) ,cellmat,cellconstr)
test_array(a,r)

a = attach_constraints_cols(cellmatvec,cellconstr)
r = map( (mat,vec,constr) -> (mat*transpose(constr),vec),cellmat,cellvec,cellconstr)
test_array(a,r)

##cache = array_cache(a)
##using BenchmarkTools
##@btime getindex!($cache,$a,2)
#
#ndofs = 4
#
#z = 2.0
#v = VectorValue(3.0,1.5)
#w = VectorValue(3.4,3.5)
#a = MockBasis{2}(v,ndofs)
#b = MockBasis{2}(w,ndofs)
#c = fill(1.0,ndofs)
#f = OtherMockBasis{2}(ndofs)
#g = MockField{2}(v)
#
#l = 10
#zl = [ z for  i in 1:l]
#cl = fill(c,l)
#fl = Fill(f,l)
#ϕl = lincomb(fl,cl)
#gl = fill(g,l)
#al = Fill(a,l)
#bl = fill(b,l)
#
#cell_axes = Fill((Base.OneTo(ndofs),),l)
#gf = GenericCellField(gl,ϕl,Val(true))
#af = GenericCellField(al,ϕl,Val(true),cell_axes,Val((:,)))
#bf = GenericCellField(bl,ϕl,Val(true),cell_axes,Val((:,)))
#zf = convert_to_cell_field(zl,ϕl)
#df = af*zf
#dft = trialize_cell_basis(df)
#
#degree = 3
#quad = CellQuadrature(degree,[QUAD,],ones(Int,l))
#
#ndofs_c = 5
#cellconstr = [ rand(ndofs_c,ndofs) for i in 1:l]
#cellvals = [ rand(ndofs) for i in 1:l]
#cellvec = integrate(bf⋅v,ϕl,quad)
#cellmat = integrate(∇(bf)⊙∇(dft),ϕl,quad)
#cellmatvec = pair_arrays(cellmat,cellvec)
#cellmatvec = attach_dirichlet(cellmatvec,cellvals)
#cellmatvec = attach_constraints_rows(cellmatvec,cellconstr)
#cellmatvec = attach_constraints_cols(cellmatvec,cellconstr)
#
#@test size(cellmatvec[1][1]) == (ndofs_c,ndofs_c)
#@test size(cellmatvec[1][2]) == (ndofs_c,)
#
#cellconstr = identity_constraints(cell_axes)
#test_array(cellconstr,Fill(Matrix(I,ndofs,ndofs),l))
#@test isa(cellconstr,Fill)
#
## Test at skeleton
#aS = merge_cell_fields_at_skeleton(af,af)
#dSt = merge_cell_fields_at_skeleton(dft,dft)
#
#axesL = Fill((Base.OneTo(ndofs),),l)
#axesR = axesL
#cellvalsL = [ rand(ndofs) for i in 1:l]
#cellvalsR = cellvalsL
#cellvals = merge_cell_dofs_at_skeleton(cellvalsL,cellvalsR,axesL,axesR)
#cellconstrL = [ rand(ndofs_c,ndofs) for i in 1:l]
#cellconstrR = cellconstrL
#axesL_rows = Fill((Base.OneTo(ndofs_c),),l)
#axesR_rows = axesL_rows
#axesL_cols = axesL
#axesR_cols = axesR
#cellconstr = merge_cell_constraints_at_skeleton(cellconstrL,cellconstrR,axesL_rows,axesR_rows,axesL_cols,axesR_cols)
#
#cellvec = integrate(jump(aS⋅v),ϕl,quad)
#cellmat = integrate( jump(aS⋅v)*(w⋅dSt.⁻),ϕl,quad)
#cellmatvec = pair_arrays(cellmat,cellvec)
#cellmatvec = attach_dirichlet(cellmatvec,cellvals)
#cellmatvec = attach_constraints_rows(cellmatvec,cellconstr)
#cellmatvec = attach_constraints_cols(cellmatvec,cellconstr)
#test_array(cellmatvec,collect(cellmatvec))
#
#cellconstr = identity_constraints(cell_axes)
#cellconstr = merge_cell_constraints_at_skeleton(cellconstr,cellconstr,cell_axes,cell_axes,cell_axes,cell_axes)
#@test isa(cellconstr,VectorOfBlockArrayCoo)
#
###a = cellmatvec
##a = cellconstr
##cache = array_cache(a)
##using BenchmarkTools
##@btime getindex!($cache,$a,2)

end # module


text:
module FactorNodeTest

using Test
using ForneyLab: FactorGraph, FactorNode, Clamp, Terminal, Variable, Interface, PointMass, GaussianMixture, Nonlinear, Unscented, Sampling, Extended, MomentConstraint, ChanceConstraint
using InteractiveUtils: subtypes

@testset "FactorNode" begin
    g = FactorGraph()

    # Collect all node types
    node_types = Type[Nonlinear{Unscented}, Nonlinear{Sampling}, Nonlinear{Extended}] # Initialze with non-concrete types
    stack = Type[FactorNode]

    while !isempty(stack)
        for node_type in subtypes(pop!(stack))
            if isconcretetype(node_type)
                (node_type == PointMass) && continue # skip PointMass
                push!(node_types, node_type)
            else
                push!(stack, node_type)
            end
        end
    end

    for node_type in node_types
        # Instantiate a test node
        if node_type == Clamp
            test_node = Clamp(Variable(), 0.0)
        elseif node_type == Terminal
            test_node = Terminal(Variable(), Clamp(Variable(), 0.0).interfaces[1])
        elseif node_type == GaussianMixture # Required for Vararg argument
            test_node = GaussianMixture(Variable(), Variable(), Variable(), Variable(), Variable(), Variable())
        elseif node_type <: Nonlinear
            test_node = node_type(Variable(), Variable(), g=()->())
        elseif node_type == MomentConstraint
            test_node = MomentConstraint(Variable(), g=()->(), G=0.0, eta_init=0.0)
        elseif node_type == ChanceConstraint
            test_node = ChanceConstraint(Variable(), G=(), epsilon=0.05)
        else
            constructor_argument_length = length(first(methods(node_type)).sig.parameters) - 1
            vars = [Variable() for v = 1:constructor_argument_length]
            test_node = node_type(vars...)
        end

        # Node fields should be of correct types
        @test isa(test_node.id, Symbol)
        @test isa(test_node.interfaces, Vector{Interface})
        @test isa(test_node.i, Dict)

        # Node constructor should automatically assign an id and check uniqueness
        @test !isempty(string(test_node.id))
        if isa(node_type, Clamp)
            @test_throws Exception Clamp(Variable(), 0.0; id=test_node.id)
        else
            @test_throws Exception node_type(vars...; id=test_node.id)
        end

        # Node constructor should assign interfaces to itself
        for iface in test_node.interfaces
            @test ===(iface.node, test_node)
        end

        # Node constructor should add node to graph
        @test ===(g.nodes[test_node.id], test_node)

    end

    # isless should be defined for nodes
    @test isless(Clamp(Variable(), 0.0; id=:a), Clamp(Variable(), 0.0; id=:b))
end

end #module


text:
using FunctionOperators, LinearAlgebra, Test
@testset "Auxiliary.jl" begin
    @testset "normalizeExpression" begin
        normE(str) = FunctionOperators.normalizeExpression(str)
        @test normE("A * (B * C * D)' + (A * (B' * C * D)')'") ==
                normE("(B' * C * D) * A' + A * (B * C * D)'") ==
                normE("B' * C * D * A' + A * D' * C' * B'")
        @test normE("A * (A + ((B -C) -E) + (e - f))") ==
                normE("A * (((B -E) -C)+A + (e - f))") ==
                normE("A * ((e - f) + A + ((B - E) - C))")
        @test normE("q + (A + ((B -C) -E) + (e - f)) * W") ==
                normE("(A + ((B -C) -E) + (e - f)) * W + q") ==
                normE("q + (e * W - f * W) + A * W + ((B * W - E * W) - C * W)")
        @test normE("A * (A + ((B -C) -E) + (e - f)) * W") ==
                normE("A * (((B -E) -C) * W + A * W + (e - f) * W)") ==
                normE("A * ((e * W - f * W) + A * W + ((B * W - E * W) - C * W))")
        @test normE("(A + B) * (C + D)") == "B * (D + C) + A * (D + C)"
        @test normE("(A + B) * (C * D)") == normE("(A + B) * C * D") ==
                normE("B * C * D + A * C * D")
        @test normE("(a + b) * (A + B) * (C + D)") ==
                normE("b * (B * (D + C) + A * (D + C)) + a * (B * (D + C) + A * (D + C))")
        @test normE("A + B + (E + D) + C") == normE("A + B + E + D + C")
        @test normE("(A - B) * C * D * E") == normE("A * C * D * E - B * C * D * E")
        @test normE("A * (D + C + B) * D * E * F") ==
                normE("A * (D * D * E * F + C * D * E * F + B * D * E * F)")
        @test normE("A * B * C * (E - F) * (D + E)") ==
                normE("A * B * C * (E * (D + E) - F * (D + E))")
        @test normE("A - (B - C)") == normE("A - B + C") == normE("C + A - B")
        # Well, I need to completely reorganize my code to support the following case, so I don't support it :(
        @test normE("A + B + (C - E)") == normE("(A + B + C) - E")
        @test normE("-(C - E)") == normE("-C + E")
    end
    @testset "Equality" begin
        Op₁ = FunctionOperator{Float64}("Op₁", x -> x, x -> x, (1,), (1,))
        Op₂ = FunctionOperator{Float64}("Op₂", x -> x, x -> x, (1,), (1,))
        @test Op₁ != Op₂
        @test Op₁ + Op₂ != Op₂
        @test Op₁ != Op₂ + Op₂
        @test Op₁ + Op₂ != Op₂ + Op₂
        @test Op₁ * Op₂ != Op₂ * Op₁
        @test Op₁ == Op₁
        @test Op₁ + Op₂ == Op₂ + Op₁
        @test Op₁ * Op₂ == Op₁ * Op₂
    end
    @testset "SetPlan" begin
        Op₁ = FunctionOperator{Float64}("Op₁", x -> x, x -> x, (1,), (1,))
        tricky = (Op₁ + 2.5I)'
        @test_throws ErrorException tricky * [5.]
        setPlan(tricky, (buffer, x) -> broadcast!(*, buffer, x, 2), "2x")
        @test tricky * [5.] == [10.]
        @test setPlan((Op₁ + 2.5I)', (buffer, x) -> broadcast!(*, buffer, x, 2)) * [5.] == [10.]
        @test repr("text/plain", tricky) == "FunctionOperatorComposite with eltype Float64\n    Name: ((Op₁ + (2.5*I)))'\n    Input dimensions: (1,)\n    Output dimensions: (1,)\n    Plan: 2x"
        setPlan(tricky, (buffer, x) -> broadcast!(*, buffer, x, 2))
        @test repr("text/plain", tricky) == "FunctionOperatorComposite with eltype Float64\n    Name: ((Op₁ + (2.5*I)))'\n    Input dimensions: (1,)\n    Output dimensions: (1,)\n    Plan: manually defined"
    end
    @testset "Macro" begin
        @testset "🔝 marker" begin
            result, var1, var2, var3 = rand(4)
            @♻ for i=1:5
                result = 🔝(var1 + var2) * var3
            end
            @test 🔝_1 == var1 + var2
            @test result == (var1 + var2) * var3
        end
        @testset "🔃 marker" begin
            result, var1, var2, var3 = rand(3,3), rand(3,3), rand(3,3), rand(3,3)
            @♻ for i=1:5
                var2 .= rand(3,3)
                result .= var1 * 🔃(var2 + var3)
                @test 🔃_1 == var2 + var3
            end
            @test result == var1 * (var2 + var3)
            @♻ for i=1:5
                result = var1 * 🔃(var2 - var3)
                @test 🔃_1 == var2 - var3
            end
            @test result == var1 * (var2 - var3)
            @♻ for i=1:5
                result = var1 * 🔃(var2 * var3)
                @test 🔃_1 == var2 * var3
            end
            @test result == var1 * (var2 * var3)
            @♻ for i=1:5
                result = var1 * 🔃(var1 + (var2 * var3))
                @test 🔃_1 == var1 + (var2 * var3)
            end
            @test result == var1 * (var1 + (var2 * var3))
        end
        @testset "@🔃 marker" begin
            var1, var2 = rand(3,3), rand(3,3)
            @♻ for i=1:5
                @🔃 result = var1 * var2
            end
            @test result == var1 * var2
            var1, var2 = rand(3,3), rand(3,3)
            @♻ for i=1:5
                @🔃 result .= var1 * var2
            end
            @test result == var1 * var2
        end
        @testset "nesting" begin
            result, var1, var2, var3 = rand(3,3), rand(3,3), rand(3,3), rand(3,3)
            @♻ for i=1:5
                @🔃 result = 🔃(🔝(var1 + var2) * var3)
                @test 🔃_1 == (var1 + var2) * var3
            end
            @test 🔝_1 == var1 + var2
            @test result == (var1 + var2) * var3
        end
        @testset "@🔃 assertion" begin
            @test_throws LoadError eval(quote
               @♻ for i=1:5
                   @🔃 result = var1 + var3
               end
           end)
        end
    end
end;


text:
using DynamicHMC, LogDensityProblems
using Distributions, Parameters, Random
using ForwardDiff, StatsPlots

include(string(pwd(),"/src/simulation.jl"))
include(string(pwd(),"/gali_bayesian.jl"))
include(string(pwd(),"/mcmc/priors_v2.jl"))
include(string(pwd(),"/mcmc/foos_v2.jl"))

yy,shocks = simulate_dsge(GAMMA_0,GAMMA_1,PSI,PI,500)

prob = DSGE_Model(yy[:,2],1/3)

prob((bet = 0.99,epsilon = 6,theta=2/3,sig=1,s2=1,phi=1,phi_pi=1.5,phi_y=0.5/4,rho_v=0.5))

t = problem_transform(prob)

P = TransformedLogDensity(t,prob)

aa = [0.99,6,2/3,1,1,1,1.5,0.5/4,0.5]

true_vals = TransformVariables.inverse(t,(bet = 0.99,epsilon = 6,theta=2/3,sig=1,s2=1,phi=1,phi_pi=1.5,phi_y=0.5/4,rho_v=0.5))

LogDensityProblems.logdensity(P,true_vals)

function LogDensityProblems.capabilities(::TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:bet, :epsilon, :theta, :sig, :s2, :phi, :phi_pi, :phi_y, :rho_v),Tuple{TransformVariables.ScaledShiftedLogistic{Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ScaledShiftedLogistic{Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ScaledShiftedLogistic{Float64}}}},DSGE_Model{Array{Float64,1},Float64}})
    LogDensityProblems.LogDensityOrder{1}() # can do gradient
end

#LogDensityProblems.dimension(::NormalPosterior) = 2 # for this problem

function LogDensityProblems.logdensity_and_gradient(problem::TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:bet, :epsilon, :theta, :sig, :s2, :phi, :phi_pi, :phi_y, :rho_v),Tuple{TransformVariables.ScaledShiftedLogistic{Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ScaledShiftedLogistic{Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ShiftedExp{true,Float64},TransformVariables.ScaledShiftedLogistic{Float64}}}},DSGE_Model{Array{Float64,1},Float64}}, x)
    logdens, grad = LogDensityProblems.logdensity(P,x)
    logdens, grad
end

res = @time mcmc_with_warmup(Random.GLOBAL_RNG, P,1000)

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC


text:
using Test
using DocumenterTools

@testset "DocumenterTools" begin
    @testset "generate" begin
        include("generate.jl")
    end

    @testset "genkeys-added" begin
        using Example
        DocumenterTools.genkeys(user="JuliaLang", repo="git@github.com:JuliaLang/Example.jl.git")
    end

    @testset "genkeys-deved" begin
        using Pkg
        Pkg.develop("DocumenterMarkdown")
        using DocumenterMarkdown
        DocumenterTools.genkeys(DocumenterMarkdown)
    end

    @testset "outdated warnings" begin
        include("outdated.jl")
    end
end


text:
# include(joinpath(dirname(dirname(@__FILE__)), "src", "SpelledOut.jl")); using .SpelledOut
using Test
using SpelledOut

@testset "English" begin
    @test spelled_out(1234, lang = :en) == "one thousand, two hundred thirty-four"
    @test spelled_out(1234, lang = :en_US) == "one thousand, two hundred thirty-four"
    @test spelled_out(1234, lang = :en_UK) == "one thousand, two hundred and thirty-four"
    @test spelled_out(123456789, lang = :en) == "one hundred twenty-three million, four hundred fifty-six thousand, seven hundred eighty-nine"
    @test spelled_out(123456789, lang = :en_UK) == "one hundred and twenty-three million, four hundred and fifty-six thousand, seven hundred and eighty-nine"
    @test spelled_out(0, lang = :en) == "zero"
    @test spelled_out(0, lang = :en_UK) == "zero"
    @test spelled_out(100, lang = :en) == "one hundred"
    @test spelled_out(100, lang = :en_UK) == "one hundred"
    @test spelled_out(123456789, lang = :en_UK, dict=:european) == "one hundred and twenty-three million, four hundred and fifty-six thousand, seven hundred and eighty-nine"
    @test spelled_out(123456789, lang = :en, dict=:european) == "one hundred twenty-three million, four hundred fifty-six thousand, seven hundred eighty-nine"
    @test spelled_out(123456789, lang = :en_UK, dict=:british) == "one hundred and twenty-three million, four hundred and fifty-six thousand, seven hundred and eighty-nine"
    @test spelled_out(1234567890123, lang = :en_UK, dict=:british) == "one billion, two hundred and thirty-four thousand million, five hundred and sixty-seven million, eight hundred and ninety thousand, one hundred and twenty-three"
    @test spelled_out(1234567890123, lang = :en_UK, dict=:modern) == "one trillion, two hundred and thirty-four billion, five hundred and sixty-seven million, eight hundred and ninety thousand, one hundred and twenty-three"
    @test spelled_out(3.0, lang = :en) == "three"
    @test spelled_out(3, lang = :en) == "three"
    @test spelled_out(3.141592653, lang = :en) == "three point one four one five nine two six five three"
    @test spelled_out(2^log2(3390000), lang = :en) == "three million, three hundred ninety thousand" # should convert 3.3899999999999986e6 to an integer
end

@testset "Spanish" begin
    @test spelled_out(585000, lang = :es) == "quinientos ochenta y cinco mil"

end

nothing


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

# The version of this JLL is decoupled from the upstream version.
# Whenever we package a new upstream release, we initially map its
# version X.Y.Z to X00.Y00.Z00 (i.e., multiply each component by 100).
# So for example version 2.6.3 would become 200.600.300.
#
# Moreover, all our packages using this JLL use `~` in their compat ranges.
#
# Together, this allows us to increment the patch level of the JLL for minor tweaks.
# If a rebuild of the JLL is needed which keeps the upstream version identical
# but breaks ABI compatibility for any reason, we can increment the minor version
# e.g. go from 200.600.300 to 200.601.300.
# To package prerelease versions, we can also adjust the minor version; e.g. we may
# map a prerelease of 2.7.0 to 200.690.000.
#
# There is currently no plan to change the major version, except when upstream itself
# changes its major version. It simply seemed sensible to apply the same transformation
# to all components.
#

# WARNING WARNING WARNING: any change to the the version of this JLL should be carefully
# coordinated with corresponding changes to Singular_jll.jl, LoadFlint.jl, Nemo.jl,
# and possibly other packages.
name = "FLINT"
upstream_version = v"2.8.4"
build_for_julia16_or_newer = true
version_offset = build_for_julia16_or_newer ? v"0.0.1" : v"0.0.0"
version = VersionNumber(upstream_version.major * 100 + version_offset.major,
                        upstream_version.minor * 100 + version_offset.minor,
                        upstream_version.patch * 100 + version_offset.patch)

# Collection of sources required to build FLINT
sources = [
    GitSource("https://github.com/wbhart/flint2.git", "432ce475efb1fa846518abf98a182ff5fa3bce4c"), # git tag v2.8.4
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/flint*
if [[ ${target} == *musl* ]]; then
   # because of some ordering issue with pthread.h and sched.h includes
   export CFLAGS=-D_GNU_SOURCE
elif [[ ${target} == *mingw* ]]; then
   extraflags=--reentrant
fi

for f in ${WORKSPACE}/srcdir/patches/*.patch; do
  atomic_patch -p1 ${f}
done

./configure --prefix=$prefix --disable-static --enable-shared --with-gmp=$prefix --with-mpfr=$prefix --with-blas=$prefix ${extraflags}
make -j${nproc}
make install LIBDIR=$(basename ${libdir})
"""

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

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

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("GMP_jll", build_for_julia16_or_newer ? v"6.2.0" : v"6.1.2"),
    Dependency("MPFR_jll", build_for_julia16_or_newer ? v"4.1.1" : v"4.0.2"),
    Dependency("OpenBLAS32_jll", v"0.3.10"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               julia_compat = build_for_julia16_or_newer ? "1.6" : "~1.0, ~1.1, ~1.2, ~1.3, ~1.4, ~1.5",
               init_block = """
  if !Sys.iswindows() && !(get(ENV, "NEMO_THREADED", "") == "1")
    #to match the global gmp ones
    fm = dlsym(libflint_handle, :__flint_set_memory_functions)
    ccall(fm, Nothing,
      (Ptr{Nothing},Ptr{Nothing},Ptr{Nothing},Ptr{Nothing}),
        cglobal(:jl_malloc),
        cglobal(:jl_calloc),
        cglobal(:jl_realloc),
        cglobal(:jl_free))
  end
""")

max_stars_repo_name: jakebolewski/Yggdrasil

max_issues_repo_name: jakebolewski/Yggdrasil

max_forks_repo_name: jakebolewski/Yggdrasil


text:
"""
```
create_q4q4_mb(mb)
```

Given a `MeansBands` object with `product in [:hist4q, :histforecast4q, :forecast4q, :bddforecast4q]`,
returns a `MeansBands` object with only Q4 observations.
"""
function create_q4q4_mb(mb::MeansBands)
    prod = mb.metadata[:product]
    if !(prod in [:hist4q, :histforecast4q, :forecast4q, :bddforecast4q])
        error("Must pass in MeansBands with product in [:hist4q, :histforecast4q, :forecast4q, :bddforecast4q] to create_forecastq4q4")
    end

    isq4(date::Date) = Dates.quarterofyear(date) == 4

    # Update metadata
    q4_metadata  = deepcopy(mb.metadata)
    q4_dates     = sort(Base.filter(isq4, collect(keys(mb.metadata[:date_inds]))))
    q4_metadata[:date_inds] = Dict(d::Date => i::Int for (i, d) in enumerate(q4_dates))
    q4_metadata[:product]   = Symbol(replace(string(prod), r"(.)4q" => s"\g<1>q4q4")) # Replace "4q" with "q4q4" in product

    # Index out Q4 observations
    q4_means = mb.means[isq4.(mb.means[!, :date]), :]

    q4_bands = Dict{Symbol, DataFrame}()
    for var in keys(mb.bands)
        dict = mb.bands[var]
        q4_bands[var] = dict[isq4.(dict[!, :date]), :]
    end

    # Return new MeansBands object
    return MeansBands(q4_metadata, q4_means, q4_bands)
end


text:
@testset "1281.subtract-the-product-and-sum-of-digits-of-an-integer.jl" begin
    @test subtract_prod_sum(234) == 15
    @test subtract_prod_sum(4421) == 21
end


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

const Mjprec =
    SMatrix{6,6}(+0.9999256782,     +0.0111820610,     +0.0048579479,
                 -0.000551,         +0.238514,         -0.435623,
                 -0.0111820611,     +0.9999374784,     -0.0000271474,
                 -0.238565,         -0.002667,         +0.012254,
                 -0.0048579477,     -0.0000271765,     +0.9999881997,
                 +0.435739,         -0.008541,         +0.002117,
                 +0.00000242395018, +0.00000002710663, +0.00000001177656,
                 +0.99994704,       +0.01118251,       +0.00485767,
                 -0.00000002710663, +0.00000242397878, -0.00000000006582,
                 -0.01118251,       +0.99995883,       -0.00002714,
                 -0.00000001177656, -0.00000000006587, 0.00000242410173,
                 -0.00485767,       -0.00002718,       1.00000956)

# Note: IDL version of `jprecess' changes in-place "muradec", "parallax" and
# "radvel".  We don't do anything like this, but calculations are below,
# commented, in case someone is interested.
function _jprecess(ra::T, dec::T, parallax::T, radvel::T, epoch::T,
                   muradec::Vector{T}) where {T<:AbstractFloat}
    @assert length(muradec) == 2
    sinra, cosra  = sincos(deg2rad(ra))
    sindec, cosdec = sincos(deg2rad(dec))
    if isfinite(epoch) && epoch != 1950
        A = A_precess .+ deg2rad.(A_dot_precess .* (epoch .- 1950) ./ 360000)
    else
        A = A_precess
    end
    r0 = SVector(cosra*cosdec,  sinra*cosdec,  sindec)
    r0_dot = SVector(-muradec[1]*sinra*cosdec - muradec[2]*cosra*sindec,
                     muradec[1]*cosra*cosdec - muradec[2]*sinra*sindec,
                     muradec[2]*cosdec) + 21.095 * radvel * parallax * r0
    r1 = r0 .- A .+ dot(r0, A) .* r0
    r1_dot = r0_dot .- A_dot_precess .+ dot(r0, A_dot_precess) .* r0
    R_1 = vcat(r1, r1_dot)
    R  = Mjprec*R_1
    if isfinite(epoch)
        t  = ((epoch - 1950) - 50.00021) / 100
        rr1 = @view(R[1:3]) .+ deg2rad.(@view(R[4:6]) .* t ./ 3600)
        x = rr1[1]
        y = rr1[2]
        z = rr1[3]
    else
        x = R[1]
        y = R[2]
        z = R[3]
        # x_dot = R[3]
        # y_dot= R[4]
        # z_dot = R[5]
    end
    rmag = norm((x, y, z))
    r2 = rmag*rmag
    dec2000 = asin(z / rmag)
    ra2000  = atan(y, x)
    # if isnan(epoch)
    #     muradec[1] = (x*y_dot - y*x_dot)/(x*x + y*y)
    #     muradec[2] = (z_dot*(x*x + y*y) - z*(x*x_dot + y*y_dot))/(r2*norm((x, y)))
    # end
    # if parallax > 0
    #     radvel   = (x*x_dot + y*y_dot + z*z_dot)/(21.095*parallax*rmag)
    #     parallax = parallax / rmag
    # end
    if ra2000 < 0
        ra2000 += 2 * T(pi)
    end
    ra2000  = rad2deg(ra2000)
    dec2000 = rad2deg(dec2000)
    return ra2000, dec2000
end

# Main interface.
jprecess(ra::Real, dec::Real, muradec::Vector{R};
         parallax::Real=0.0, radvel::Real=0.0) where {R<:Real} =
             _jprecess(promote(float(ra), float(dec), float(parallax),
                               float(radvel), NaN)...,
                       float(muradec))

jprecess(ra::Real, dec::Real, epoch::Real=1950.0) =
    _jprecess(promote(float(ra), float(dec), 0.0, 0.0, float(epoch))...,
              zeros(typeof(float(ra)), 2))

# Tuple arguments.
jprecess(radec::Tuple{Real,Real}, muradec::Vector{<:Real};
         parallax::Real=0.0, radvel::Real=0.0) =
             jprecess(radec..., muradec, parallax=parallax, radvel=radvel)

jprecess(radec::Tuple{Real,Real}, epoch::Real=1950.0) =
    jprecess(radec..., epoch)

# Vectorial arguments.
function jprecess(ra::AbstractArray{R}, dec::AbstractArray{<:Real},
                  muradec::AbstractArray{<:Real};
                  parallax::AbstractArray{<:Real}=zeros(R, length(ra)),
                  radvel::AbstractArray{<:Real}=zeros(R, length(ra))) where {R<:Real}
    @assert length(ra) == length(dec) == size(muradec)[2] == length(parallax) == length(radvel)
    typer = float(R)
    ra2000  = similar(ra, typer)
    dec2000 = similar(dec, typer)
    for i in eachindex(ra)
        ra2000[i], dec2000[i] = jprecess(ra[i], dec[i], muradec[:,i],
                                         parallax=parallax[i], radvel=radvel[i])
    end
    return ra2000, dec2000
end

function jprecess(ra::AbstractArray{R}, dec::AbstractArray{<:Real},
                  epoch::Real=1950.0) where {R<:Real}
    @assert length(ra) == length(dec)
    typer = float(R)
    ra2000  = similar(ra, typer)
    dec2000 = similar(dec, typer)
    for i in eachindex(ra)
        ra2000[i], dec2000[i] = jprecess(ra[i], dec[i], epoch[i])
    end
    return ra2000, dec2000
end

"""
    jprecess(ra, dec[, epoch]) -> ra2000, dec2000
    jprecess(ra, dec, muradec[, parallax=parallax, radvel=radvel]) -> ra2000, dec2000

### Purpose ###

Precess positions from B1950.0 (FK4) to J2000.0 (FK5).

### Explanation ###

Calculate the mean place of a star at J2000.0 on the FK5 system from the mean
place at B1950.0 on the FK4 system.

`jprecess` function has two methods, one for each of the following cases:

* the proper motion is known and non-zero
* the proper motion is unknown or known to be exactly zero (i.e. extragalactic
  radio sources).  Better precision can be achieved in this case by inputting
  the epoch of the original observations.

### Arguments ###

The function has 2 methods.  The common mandatory arguments are:

* `ra`: input B1950 right ascension, in degrees.
* `dec`: input B1950 declination, in degrees.

The two methods have a different third argument (see "Explanation" section for
more details).  It can be one of the following:

* `muradec`: 2-element vector containing the proper motion in seconds of arc
  per tropical *century* in right ascension and declination.
* `epoch`: scalar giving epoch of original observations.

If none of these two arguments is provided (so `jprecess` is fed only with right
ascension and declination), it is assumed that proper motion is exactly zero and
`epoch = 1950`.

If it is used the method involving `muradec` argument, the following keywords
are available:

* `parallax` (optional numerical keyword): stellar parallax, in seconds of arc.
* `radvel` (optional numerical keyword): radial velocity in km/s.

Right ascension and declination can be passed as the 2-tuple `(ra, dec)`.  You
can also pass `ra`, `dec`, `parallax`, and `radvel` as arrays, all of the same
length N.  In that case, `muradec` should be a matrix 2×N.

### Output ###

The 2-tuple of right ascension and declination in 2000, in degrees, of input
coordinates is returned.  If `ra` and `dec` (and other possible optional
arguments) are arrays, the 2-tuple of arrays `(ra2000, dec2000)` of the same
length as the input coordinates is returned.

### Method ###

The algorithm is taken from the Explanatory Supplement to the Astronomical
Almanac 1992, page 184.  See also Aoki et al (1983), A&A, 128, 263.  URL:
http://adsabs.harvard.edu/abs/1983A%26A...128..263A.

### Example ###

The SAO catalogue gives the B1950 position and proper motion for the
star HD 119288.   Find the J2000 position.

* RA(1950) = 13h 39m 44.526s
* Dec(1950) = 8d 38' 28.63''
* Mu(RA) = -.0259 s/yr
* Mu(Dec) = -.093 ''/yr

```jldoctest
julia> using AstroLib

julia> muradec = 100 * [-15*0.0259, -0.093]; # convert to century proper motion

julia> ra = ten(13, 39, 44.526)*15;

julia> decl = ten(8, 38, 28.63);

julia> adstring(jprecess(ra, decl, muradec), precision=2)
" 13 42 12.740  +08 23 17.69"
```

### Notes ###

"When transferring individual observations, as opposed to catalog mean place,
the safest method is to tranform the observations back to the epoch of the
observation, on the FK4 system (or in the system that was used to to produce the
observed mean place), convert to the FK5 system, and transform to the the epoch
and equinox of J2000.0" -- from the Explanatory Supplement (1992), p. 180

`bprecess` performs the precession to B1950 coordinates.

Code of this function is based on IDL Astronomy User's Library.
"""
jprecess


text:
using WignerFamilies
using HalfIntegers
# using OffsetArrays

# wsv = WignerSymbolVector(zeros(length(HalfInt(3/2):HalfInt(7/2))), HalfInt(3/2):HalfInt(7/2))

w = WignerF( HalfInt(5/2), 5, HalfInt(1/2), HalfInt(-1) )
m₁ = -w.m₂ - w.m₃

w3j = get_wigner_array(w)
wigner3j_f!(w, w3j)

using WignerSymbols
# @test wigner3j(Float64, HalfInt(5/2), w.j₂, w.j₃, m₁, w.m₂ ) ==

##


text:
using Test
using SafeTestsets

@testset "advection" begin
  @safetestset "wave_1d" begin include("wave_1d.jl") end
  @safetestset "wave_2d" begin include("wave_2d.jl") end
  @safetestset "wave_3d" begin include("wave_3d.jl") end
  @safetestset "entropy_conservation_1d" begin include("entropy_conservation_1d.jl") end
end


text:
# This file is a part of ArraysOfArrays.jl, licensed under the MIT License (MIT).

@inline flatview(A::AbstractArray{SA,N}) where {S,T,M,N,SA<:StaticArrays.StaticArray{S,T,M}} =
    reshape(reinterpret(T, A), size(SA)..., size(A)...)

@inline function nestedview(A::AbstractArray{T}, SA::Type{StaticArrays.SVector{S,T}}) where {T,S}
    size_A = size(A)
    size_A[1] == S || throw(DimensionMismatch("Length $S of static vector type does not match first dimension of array of size $size_A"))
    reshape(reinterpret(SA, A), _tail(size_A)...)
end

@inline nestedview(A::AbstractArray{T}, ::Type{StaticArrays.SVector{S}}) where {T,S} =
    nestedview(A, StaticArrays.SVector{S,T})


text:
using DataFrames
using HDF5

function main()
    table_path = ARGS[1]
    id = h5read(table_path,"id")
    x = h5read(table_path,"x")
    y = h5read(table_path,"y")

    df = DataFrame()
    df[:id] = id
    df[:x] = x
    df[:y] = y
    tic()
    df[df[:id].<100,:]
    t = toq()
    println("DataFrames.jl filter took: ", t)
end

main()

max_stars_repo_name: ehsantn/data-analytics-benchmarks

max_issues_repo_name: ehsantn/data-analytics-benchmarks

max_forks_repo_name: ehsantn/data-analytics-benchmarks


text:
abstract type AbstractVirtualVolume{T} end

in(p::AbstractCoordinatePoint{T, 3}, avv::AbstractVirtualVolume{T}) where {T <: SSDFloat} = in(p, avv.geometry)


text:
@testset "MInt" begin

@test Metaprogramming.MInt_from_VoidTuple(()) === MInt{0}
@test Metaprogramming.MInt_from_VoidTuple((nothing,)) === MInt{1}
@test Metaprogramming.MInt_from_VoidTuple((nothing,nothing,nothing)) === MInt{3}

@test Metaprogramming.VoidTuple_from_MInt(MInt{0}) === ()
@test Metaprogramming.VoidTuple_from_MInt(MInt{1}) === (nothing,)
@test Metaprogramming.VoidTuple_from_MInt(MInt{3}) === (nothing,nothing,nothing)

@test begin
    a = MInt{2}
    b = MInt{10}
    a+b === MInt{12}
end

@test begin
    a = MInt{0}
    b = MInt{0}
    a+b === MInt{0}
end

@test begin
    a = MInt{0}
    b = MInt{3}
    a+b === MInt{3}
end

@test begin
    a = MInt{10}
    b = MInt{2}
    a-b === MInt{8}
end

@test begin
    a = MInt{0}
    b = MInt{0}
    a-b === MInt{0}
end

@test begin
    a = MInt{4}
    b = MInt{0}
    a-b === MInt{4}
end

@test begin
    a = MInt{4}
    b = MInt{4}
    a-b === MInt{0}
end

end


text:
module HCIDatasets

using DataDeps
using FITSIO

# helper for loading data with FITSIO.jl
getdata(path) = read(FITS(path)[1])

# setup type system
abstract type HCIDataset end

Base.getindex(d::HCIDataset, keys...) = map(k -> d[k], keys)
Base.getindex(d::HCIDataset, keys) = map(k -> d[k], keys)
Base.pairs(d::HCIDataset) = (k => d[k] for k in keys(d))

include("HR8799.jl")
include("BetaPictoris.jl")
include("V471Tau.jl")

function __init__()
    register(HR8799_datadep)
    register(BetaPictoris_datadep)
    register(V471Tau_datadep)
end

end


text:
using ArgParse;
using Distributions;
using LinearAlgebra;
using Logging;
using DelimitedFiles;
#using ProfileView;
using DecFP;
using Quadmath;

include(joinpath("inc", "eap_chain.jl"));
include(joinpath("inc", "average.jl"));
include(joinpath("inc", "acceptance.jl"));

s = ArgParseSettings();
@add_arg_table! s begin
  "--E0", "-e"
    help = "magnitude of electric field"
    arg_type = Float64
    default = 0.0
  "--chain-type", "-T"
    help = "chain type (dielectric|polar)"
    arg_type = String
    default = "dielectric"
  "--K1", "-J"
    help = "dipole susceptibility along the monomer axis (dielectric chain)"
    arg_type = Float64
    default = 1.0
  "--K2", "-K"
    help = "dipole susceptibility orthogonal to the monomer axis (dielectric chain)"
    arg_type = Float64
    default = 0.0
  "--mu", "-m"
    arg_type = Float64
    default = 1e-2
    help = "dipole magnitude (electret chain)"
  "--energy-type", "-u"
    help = "energy type (noninteracting|interacting|Ising)"
    arg_type = String
    default = "noninteracting"
  "--kT", "-k"
    help = "dimensionless temperature"
    arg_type = Float64
    default = 1.0
  "--Fz", "-F"
    help = "force in the z-direction (direction of E-field; force ensemble)"
    arg_type = Float64
    default = 0.0
  "--Fx", "-G"
    help = "force in the x-direction (force ensemble)"
    arg_type = Float64
    default = 0.0
  "--mlen", "-b"
    help = "monomer length"
    arg_type = Float64
    default = 1.0
  "--num-monomers", "-n"
    help = "number of monomers"
    arg_type = Int
    default = 100
  "--num-steps", "-N"
    help = "number of steps"
    arg_type = Int
    default = convert(Int, 1e6)
  "--phi-step", "-p"
    help = "maximum ϕ step length"
    arg_type = Float64;
    default = 3*π / 8;
  "--theta-step", "-q"
    help = "maximum θ step length"
    arg_type = Float64;
    default = 3*π / 16;
#  "--do-clustering"
#    help = "trial moves with clustering"
#    action = :store_true
  "--cluster-prob"
    help = "probability of flipping a cluster"
    arg_type = Float64;
    default = 0.5;
  "--step-adjust-lb", "-L"
    help = "adjust step sizes if acc. ratio below this threshold"
    arg_type = Float64
    default = 0.15
  "--step-adjust-ub", "-U"
    help = "adjust step sizes if acc. ratio above this threshold"
    arg_type = Float64
    default = 0.40
  "--step-adjust-scale", "-A"
    help = "scale factor for adjusting step sizes (> 1.0)"
    arg_type = Float64
    default = 1.1
  "--steps-per-adjust", "-S"
    help = "steps between step size adjustments"
    arg_type = Int
    default = 2500
  "--umbrella-sampling", "-B"
    help = "use umbrella sampling (w/ electrostatic weight function)"
    action = :store_true
  "--update-freq"
    help = "update frequency (seconds)"
    arg_type = Float64;
    default = 15.0;
  "--verbose", "-v"
    help = "verbosity level: 0-nothing, 1-errors, 2-warnings, 3-info"
    arg_type = Int
    default = 3
  "--prefix", "-P"
    help = "prefix for output files"
    arg_type = String
    default = "eap-mcmc"
  "--postfix", "-Q"
    help = "postfix for output files"
    arg_type = String
    default = ""
  "--stepout", "-s"
    help = "steps between storing microstates"
    arg_type = Int
    default = 500
  "--numeric-type"
    help = "numerical data type for averaging (float64|float128|dec128|big)"
    arg_type = String
    default = "float64"
  "--profile", "-Z"
    help = "profile the program"
    action = :store_true
end

pargs = parse_args(s);

if pargs["verbose"] == 3
  global_logger(ConsoleLogger(stderr, Logging.Info));
elseif pargs["verbose"] == 2
  global_logger(ConsoleLogger(stderr, Logging.Warn));
elseif pargs["verbose"] == 1
  global_logger(ConsoleLogger(stderr, Logging.Error));
else
  global_logger(Logging.NullLogger());
end

function mcmc(nsteps::Int, pargs)
  ϕstep, θstep = pargs["phi-step"], pargs["theta-step"];
  dϕ_dist = Uniform(-ϕstep, ϕstep);
  dθ_dist = Uniform(-θstep, θstep);
  chain = EAPChain(pargs);
  chain.U = U(chain);
  wf = AntiDipoleWeightFunction(chain);
  acceptor = Metropolis(
               chain,
               (pargs["umbrella-sampling"]) ? wf : WeightlessFunction()
              );
  numeric_type = if pargs["numeric-type"] == "float64"
    Float64;
  elseif pargs["numeric-type"] == "float128"
    Float128;
  elseif pargs["numeric-type"] == "dec128"
    Dec128;
  elseif pargs["numeric-type"] == "big"
    BigFloat;
  else
    error("numeric-type '$(pargs["numeric-type"])' not understood");
    exit(-1);
  end

  # setup averagers
  Avg, avgcons, avgconss = if pargs["umbrella-sampling"]
    (UmbrellaAverager,
      ((acc, chain) -> begin;
        UmbrellaAverager(
                         StandardAverager(
                                          numeric_type(0.0*acc(chain)),
                                          numeric_type(0.0),
                                          acc
                                         ),
                         wf
                        );
     end),
     ((acc, chain) -> begin;
        UmbrellaAverager(
                         StandardAverager(
                                          Vector{numeric_type}(0.0*acc(chain)),
                                          numeric_type(0.0),
                                          acc
                                         ),
                         wf
                        );
     end)
     );
  else
    (StandardAverager,
      ((acc, chain) -> begin;
             StandardAverager(
                              numeric_type,
                              acc, chain
                             )
     end),
     ((acc, chain) -> begin;
             StandardAverager(
                              numeric_type,
                              Vector{numeric_type},
                              acc, chain
                             );
     end)
     );
  end
  scalar_averagers = (Avg{numeric_type,numeric_type}[
                       (avgcons(chain -> dot(end_to_end(chain),
                                             end_to_end(chain)), chain)),
                       (avgcons(chain -> dot(chain_μ(chain), chain_μ(chain)),
                                chain)),
                       (avgcons(chain -> chain.U, chain)),
                       (avgcons(chain -> chain.U*chain.U, chain)),
                      ]);
  vector_averagers = (Avg{Vector{numeric_type},numeric_type}[
                       avgconss(end_to_end, chain),
                       avgconss(chain -> map(x -> x*x, end_to_end(chain)), chain),
                       avgconss(chain_μ, chain),
                       avgconss(chain -> map(x -> x*x, chain_μ(chain)), chain)
                      ]);

  outfile = open("$(pargs["prefix"])_trajectory.csv", "w");
  writedlm(outfile, ["step" "r1" "r2" "r3" "p1" "p2" "p3" "U"], ',');
  rollfile = open("$(pargs["prefix"])_rolling.csv", "w");
  writedlm(rollfile, ["step" "r1" "r2" "r3" "r1sq" "r2sq" "r3sq" "rsq" "p1" "p2" "p3" "p1sq" "p2sq" "p3sq" "psq" "U" "Usq"], ',');

  start = time();
  last_update = start;
  nacc = 0;
  nacc_total = 0;
  natt = 0;

  for step=1:nsteps
    idx = rand(1:n(chain));
    dϕ = rand(dϕ_dist);
    dθ = rand(dθ_dist);
    trial_chain = EAPChain(chain);
    successful = move!(trial_chain, idx, dϕ, dθ);
    α = cluster_flip!(trial_chain, idx; ϵflip = pargs["cluster-prob"]);
    if successful && acceptor(trial_chain, rand(); α = α)
      chain = trial_chain;
      nacc += 1;
      nacc_total += 1;
    end
    natt += 1;

    if time() - last_update > pargs["update-freq"]
      @info "elapsed: $(time() - start)";
      @info "step:    $step / $nsteps";
      last_update = time();
    end

    if (
        pargs["step-adjust-scale"] != 1.0 &&
        step % pargs["steps-per-adjust"] == 0
       ) # adjust step size?
      ar = nacc / natt;
      if (ar > pargs["step-adjust-ub"] &&
          ϕstep != π && θstep != π/2)
        @info "acceptance ratio is high; increasing step size";
        nacc = 0;
        natt = 0;
        ϕstep = min(π, ϕstep*pargs["step-adjust-scale"]);
        θstep = min(π/2, θstep*pargs["step-adjust-scale"]);
      elseif ar < pargs["step-adjust-lb"]
        @info "acceptance ratio is low; decreasing step size";
        nacc = 0;
        natt = 0;
        ϕstep /= pargs["step-adjust-scale"];
        θstep /= pargs["step-adjust-scale"];
      end
      dϕ_dist = Uniform(-ϕstep, ϕstep);
      dθ_dist = Uniform(-θstep, θstep);
    end

    foreach(a -> record!(a, chain), scalar_averagers);
    foreach(a -> record!(a, chain), vector_averagers);
    if step % pargs["stepout"] == 0
      writedlm(outfile,
               hcat(step, transpose(chain.r),
                    transpose(chain_μ(chain)), chain.U),
               ',');
      writedlm(rollfile,
               hcat(
                    step,
                    transpose(get_avg(vector_averagers[1])), # r
                    transpose(get_avg(vector_averagers[2])), # rj2
                    get_avg(scalar_averagers[1]), # r2
                    transpose(get_avg(vector_averagers[3])), # p
                    transpose(get_avg(vector_averagers[4])), # pj2
                    get_avg(scalar_averagers[2]), # p2
                    get_avg(scalar_averagers[3]), # U
                    get_avg(scalar_averagers[4])  # U2
                   ),
               ',');

    end

  end # steps

  AR = nacc_total / pargs["num-steps"];
  @info "total time elapsed: $(time() - start)";
  @info "acceptance rate: $AR";

  close(outfile);
  close(rollfile);

  return (scalar_averagers, vector_averagers, AR);
end

(sas, vas, ar) = if pargs["profile"]
  @info "Profiling the mcmc code...";
  mcmc(5, pargs); # run first to compile code
  @profview mcmc(pargs["num-steps"], pargs);
  println("Press ENTER to continue...");
  readline(stdin);
  exit(1);
else
  mcmc(pargs["num-steps"], pargs);
end

println("<r>    =   $(get_avg(vas[1]))");
println("<r/nb> =   $(get_avg(vas[1]) / (pargs["mlen"]*pargs["num-monomers"]))");
println("<rj2>  =   $(get_avg(vas[2]))");
println("<r2>   =   $(get_avg(sas[1]))");
println("<p>    =   $(get_avg(vas[3]))");
println("<pj2>  =   $(get_avg(vas[4]))");
println("<p2>   =   $(get_avg(sas[2]))");
println("<U>    =   $(get_avg(sas[3]))");
println("<U2>   =   $(get_avg(sas[4]))");
println("AR     =   $ar");

max_stars_repo_name: grasingerm/polymer-stats

max_issues_repo_name: grasingerm/polymer-stats

max_forks_repo_name: grasingerm/polymer-stats