text:
function loopwinAnnot(nReads, rawdata, grpInfo, offsetSize, snpID,
                      bin2geno, nPer, SNPSize, flagAnnotate,
                      snpWtType, keepIdx, nPerKeep, kernel,
                      kinship, test, XtNullBasis, WPreSim,
                      nBlockAscent, nMMmax, nNullSimPts, tolX,
                      pvalueComputing, device, nNullSimNewtonIter,
                      rankQPhi, evalPhiAdj, XPhitNullBasis, yShift,
                      y, KPhiAdj, weightedW, QPhi, PrePartialSumW,
                      PreTotalSumW, windowSize, nPreRank, X, snpPos,
                      chrID, geneName, ynew)

  QRes = Array(Float64, nPerKeep, rankQPhi);
  tmpvec = similar(yShift);
  tmpvecQRes = Array(Float64, rankQPhi);
  yWork = similar(evalPhiAdj);
  partialSumWConst = Array(Float64, nNullSimPts);
  totalSumWConst = Array(Float64, nNullSimPts);
  subXPhitSV = Array(Float64, size(XPhitNullBasis, 2), rankQPhi);
  pSubXPhitSV = pointer(subXPhitSV);
  offset = 0;
  oneConst = ones(nPerKeep);
  geno = Array(Float64, nPer, windowSize);
  genoOri = Array(Float64, nPer, windowSize);
  snpIDred = Array(String, windowSize);
  snpPosred = Array(String, windowSize);
  mafOri = Array(Float64, windowSize);
  maf = Array(Float64, windowSize);
  sumIsnan = Array(Float64, 1, windowSize);
  nChrObs = Array(Float64, 1, windowSize);
  if pvalueComputing == "chi2"
    nPtsChi2 = 300;
    partialSumW = Array(Float64, nPtsChi2);
    totalSumW = Array(Float64, nPtsChi2);
    lambda = Array(Float64, 1, nPtsChi2);
    W = Array(Float64, windowSize, nPtsChi2);
    tmpmat0 = Array(Float64, windowSize, nPtsChi2);
    tmpmat1 = Array(Float64, windowSize, nPtsChi2);
    tmpmat2 = Array(Float64, windowSize, nPtsChi2);
    tmpmat3 = Array(Float64, windowSize, nPtsChi2);
    tmpmat4 = Array(Float64, windowSize, nPtsChi2);
    tmpmat5 = Array(Float64, windowSize, nPtsChi2);
    denomvec = Array(Float64, nPtsChi2);
    d1f = Array(Float64, nPtsChi2);
    d2f = Array(Float64, nPtsChi2);
    if test == "eScore"
      simnull = Array(Float64, nPtsChi2);
    else
      simnull = Float64[];
    end
  else
    nPtsChi2 = 300;
    partialSumW = Array(Float64, nNullSimPts);
    totalSumW = Array(Float64, nNullSimPts);
    lambda = Array(Float64, 1, nNullSimPts);
    W = Array(Float64, windowSize, nNullSimPts);
    tmpmat0 = Array(Float64, windowSize, nNullSimPts);
    tmpmat1 = Array(Float64, windowSize, nNullSimPts);
    tmpmat2 = Array(Float64, windowSize, nNullSimPts);
    tmpmat3 = Array(Float64, windowSize, nNullSimPts);
    tmpmat4 = Array(Float64, windowSize, nNullSimPts);
    tmpmat5 = Array(Float64, windowSize, nNullSimPts);
    denomvec = Array(Float64, nNullSimPts);
    d1f = Array(Float64, nNullSimPts);
    d2f = Array(Float64, nNullSimPts);
    if test == "eScore"
      simnull = Array(Float64, nNullSimPts);
    else
      simnull = Float64[];
    end
  end

  results = Array(Any, nReads, 8);
  idx = 0;

  for idxRead = 1 : nReads

    # read the binary data
    curSNPSize = grpInfo[idxRead];
    readgeno!(genoOri, curSNPSize, nPer, SNPSize, idxRead, bin2geno, rawdata,
              offset, flagAnnotate, offsetSize = offsetSize);

    tmpMatIsnan = isnan(genoOri);
    sum!(sumIsnan, tmpMatIsnan);
    for i = 1 : curSNPSize
      nChrObs[1, i] = 2 * (nPer - sumIsnan[1, i]);
    end

    # minor allele frequencies and get rid of mono-allelic SNPs
    counter = 0;
    offset = offsetSize[idxRead];
    for i = 1 : curSNPSize
      mafOri[i] = 0.0;
      for j = 1 : nPer
        if !tmpMatIsnan[j, i]
          mafOri[i] += genoOri[j, i];
        end
      end
      mafOri[i] = mafOri[i] / nChrObs[1, i];
      if mafOri[i] != 0
        counter += 1;
        maf[counter] = mafOri[i];
        snpIDred[counter] = snpID[i + offset];
        snpPosred[counter] = snpPos[i + offset];
        pGenoOri = pointer(genoOri) + (i - 1) * nPer * sizeof(Float64);
        pGeno = pointer(geno) + (counter - 1) * nPer * sizeof(Float64);
        BLAS.blascopy!(nPer, pGenoOri, 1, pGeno, 1);
      end
    end

    if counter == 0
      continue;
    end

    # impute missing genotype by expected minor allele counts
    for j = 1 : counter
      for i = 1 : nPer
        if isnan(geno[i, j])
          geno[i, j] = 2 * maf[j];
        end
      end
    end

    # get SNP weights
    grpSize = counter;
    if isempty(snpWtType)
      # constant weights
      snpSqrtWts = ones(grpSize);
    elseif snpWtType == "beta"
      # weights by beta density evaluated at MAF
      snpSqrtWts = sqrt(Distributions.pdf(Distributions.Beta(1, 25), maf[1 : grpSize]));
    elseif snpWtType == "invvar"
      # weights by inverse variance maf*(1-maf)
      snpSqrtWts = 1 ./ sqrt( maf[1 : grpSize] .*
                             (1 - maf[1 : grpSize]) );
    end

    # retrieve group genotypes
    if kernel == "GRM"
      gMAF = maf[1 : grpSize];
      gSNP = geno[keepIdx, 1 : grpSize];
      # center and scale genotype matrix by MAF
      BLAS.ger!(-2.0, oneConst, gMAF, gSNP);
      tmpvecKernel = zeros(grpSize);
      for i = 1 : grpSize
        tmpvecKernel[i] = snpSqrtWts[i] / sqrt(2 * gMAF[i] * (1 - gMAF[i]));
      end
      scale!(gSNP, tmpvecKernel);
      #gSNP = gSNP .* (snpSqrtWts ./ sqrt(2 * gMAF .* (1 - gMAF)))';
    elseif kernel == "IBS1"
      # Table 3 in notes
      gSNP = hcat(geno[keepIdx, 1 : grpSize],
                  2 - geno[keepIdx, 1 : grpSize]);
      scale!(gSNP, [snpSqrtWts; snpSqrtWts]);
    elseif kernel == "IBS2"
      # Table 2 in notes
      gSNP = hcat(geno[keepIdx, 1 : grpSize] .>= 1,
                  geno[keepIdx, 1 : grpSize] .<= 1);
      scale!(gSNP, [snpSqrtWts; snpSqrtWts]);
    elseif kernel == "IBS3"
      # Table 1 in notes
      gSNP = hcat(-2 * (geno[keepIdx, 1 : grpSize] .<= 1) + 1,
                  2 * (geno[keepIdx, 1 : grpSize] .>= 1) - 1,
                  fill(sqrt(2), nPerKeep, grpSize));
      scale!(gSNP, [snpSqrtWts; snpSqrtWts; snpSqrtWts]);
    end

    ## testing

    if kinship == "none"

      if test == "eRLRT"
        (b, vc0List, vc1List, pvalList) =
          vctest(ynew, [], XtNullBasis * gSNP, WPreSim = WPreSim, tests = test,
                 Vform = "half", nBlockAscent = nBlockAscent, nMMmax = nMMmax,
                 nNullSimPts = nNullSimPts, pvalueComputings = pvalueComputing,
                 nNullSimNewtonIter = nNullSimNewtonIter, tolX = tolX,
                 devices = device, PrePartialSumW = PrePartialSumW,
                 PreTotalSumW = PreTotalSumW, partialSumWConst = partialSumWConst,
                 totalSumWConst = totalSumWConst, windowSize = windowSize,
                 partialSumW = partialSumW, totalSumW = totalSumW,
                 lambda = lambda, W = W, nPreRank = nPreRank,
                 tmpmat0 = tmpmat0, tmpmat1 = tmpmat1, tmpmat2 = tmpmat2,
                 tmpmat3 = tmpmat3, tmpmat4 = tmpmat4, tmpmat5 = tmpmat5,
                 denomvec = denomvec, d1f = d1f, d2f = d2f, offset = offset,
                 nPtsChi2 = nPtsChi2, simnull = simnull);
      else
        (b, vc0List, vc1List, pvalList) =
          vctest(y, X, gSNP, WPreSim = WPreSim, tests = test,
                 Vform = "half", nBlockAscent = nBlockAscent, nMMmax = nMMmax,
                 nNullSimPts = nNullSimPts, pvalueComputings = pvalueComputing,
                 nNullSimNewtonIter = nNullSimNewtonIter, tolX = tolX,
                 devices = device, PrePartialSumW = PrePartialSumW,
                 PreTotalSumW = PreTotalSumW, partialSumWConst = partialSumWConst,
                 totalSumWConst = totalSumWConst, windowSize = windowSize,
                 partialSumW = partialSumW, totalSumW = totalSumW,
                 lambda = lambda, W = W, nPreRank = nPreRank,
                 tmpmat0 = tmpmat0, tmpmat1 = tmpmat1, tmpmat2 = tmpmat2,
                 tmpmat3 = tmpmat3, tmpmat4 = tmpmat4, tmpmat5 = tmpmat5,
                 denomvec = denomvec, d1f = d1f, d2f = d2f, offset = offset,
                 nPtsChi2 = nPtsChi2, simnull = simnull);
      end

    else

      # obtain some orthonormal vectors of space R^n - [QX,QPHI,S]
      # See Golub and Van Loan (1996) Algorithm 12.4.2 on p602
      tmpMat = Array(Float64, grpSize, size(XPhitNullBasis, 2));
      BLAS.gemm!('T', 'N', 1.0, gSNP, XPhitNullBasis, 0.0, tmpMat);
      (UXPhitS, svalXPhitS, VXPhitS) = svd(tmpMat, thin = false);
      rankXPhitS =
        countnz(svalXPhitS .>
                max(size(XPhitNullBasis, 2),
                    size(gSNP, 2) * eps(svalXPhitS[1])));
      pXPhitSV = pointer(VXPhitS) +
        rankXPhitS * size(XPhitNullBasis, 2) * sizeof(Float64);
      BLAS.blascopy!(size(XPhitNullBasis, 2) * rankQPhi,
                     pXPhitSV, 1, pSubXPhitSV, 1);
      BLAS.gemm!('N', 'N', 1.0, XPhitNullBasis, subXPhitSV, 0.0, QRes);

      # working Y and working non-trivial variance component
      #yWork = (PhiAdjsvd[:U]' * (yShift + KPhiAdj * (QRes' * y))) ./ sqrt(evalPhiAdj);
      BLAS.blascopy!(length(yShift), yShift, 1, tmpvec, 1);
      BLAS.gemv!('T', 1.0, QRes, y, 0.0, tmpvecQRes);
      BLAS.gemv!('N', 1.0, KPhiAdj, tmpvecQRes, 1.0, tmpvec);
      #yWork = BLAS.gemv('T', 1.0, PhiAdjsvd[:U], tmpvec) ./ sqrt(evalPhiAdj);
      BLAS.gemv!('T', 1.0, weightedW, tmpvec, 0.0, yWork);
      #for i = 1 : length(yWork)
      #  yWork[i] = yWork[i] * InvSqrtevalPhiAdj[i];
      #end
      #VWorkSqrt = PhiAdjsvd[:U]' * (QPhi' * gSNP) ./ sqrt(evalPhiAdj);
      VWorkSqrt = Array(Float64, length(evalPhiAdj), size(gSNP, 2));
      VWorkSqrt2 = Array(Float64, length(evalPhiAdj), size(gSNP, 2));
      BLAS.gemm!('T', 'N', 1.0, QPhi, gSNP, 0.0, VWorkSqrt2);
      BLAS.gemm!('T', 'N', 1.0, weightedW, VWorkSqrt2, 0.0, VWorkSqrt);
      #scale!(InvSqrtevalPhiAdj, VWorkSqrt);

      # obtain LRT and p-value
      (b, vc0List, vc1List, pvalList) =
        vctest(yWork, [], VWorkSqrt, WPreSim = WPreSim, tests = test,
               Vform = "half", nBlockAscent = nBlockAscent, nMMmax = nMMmax,
               nNullSimPts = nNullSimPts, pvalueComputings = pvalueComputing,
               nNullSimNewtonIter = nNullSimNewtonIter, tolX = tolX,
               devices = device, PrePartialSumW = PrePartialSumW,
               PreTotalSumW = PreTotalSumW, partialSumWConst = partialSumWConst,
               totalSumWConst = totalSumWConst, windowSize = windowSize,
               partialSumW = partialSumW, totalSumW = totalSumW,
               lambda = lambda, W = W, nPreRank = nPreRank,
               tmpmat0 = tmpmat0, tmpmat1 = tmpmat1, tmpmat2 = tmpmat2,
               tmpmat3 = tmpmat3, tmpmat4 = tmpmat4, tmpmat5 = tmpmat5,
               denomvec = denomvec, d1f = d1f, d2f = d2f, offset = offset,
               nPtsChi2 = nPtsChi2, simnull = simnull);

    end

    idx += 1;
    results[idx, :] = hcat(snpIDred[1], snpIDred[counter], chrID[offset + 1],
                           snpPosred[1], snpPosred[counter], geneName[idxRead],
                           curSNPSize, pvalList);

  end

  return results[1:idx, :];

end


text:
using JuLIP
using JuLIP.Potentials
using JuLIP.Testing

pairpotentials = [
   LennardJones(1.0,1.0);
   Morse(4.0,1.0,1.0);
   SWCutoff(1.0, 3.0) * LennardJones(1.0,1.0);
   SplineCutoff(2.0, 3.0) * LennardJones(1.0,1.0);
   LennardJones(1.0, 1.0) * C2Shift(2.0);
]

println("============================================")
println("  Testing pair potential implementations ")
println("============================================")
r = linspace(0.8, 4.0, 100) |> collect
push!(r, 2.0-1e-12)
for pp in pairpotentials
   println("--------------------------------")
   println(pp)
   println("--------------------------------")
   @test fdtest(pp, r, verbose=verbose)
end

print("testing shift-cutoffs: ")
V = @analytic r -> exp(r)
Vhs = V * HS(1.0)
r1 = linspace(0.0, 1.0-1e-14, 20)
r2 = linspace(1.0+1e-14, 3.0, 20)
print("HS")
@test Vhs.(r1) == exp.(r1)
@test norm(Vhs.(r2)) == 0.0
print(", V0")
V0 = V * C0Shift(1.0)
@test V0.(r1) ≈ exp.(r1) - exp(1.0)
@test norm(V0.(r2)) == 0.0
print(", V1")
V1 = V * C1Shift(1.0)
@test V1.(r1) ≈ exp.(r1) - exp(1.0) - exp(1.0) * (r1-1.0)
@test norm(V1.(r2)) == 0.0
println(", V2")
V2 = V * C2Shift(1.0)
@test V2.(r1) ≈ exp.(r1) - exp(1.0) - exp(1.0) * (r1-1.0) - 0.5 * exp(1.0) * (r1-1.0).^2
@test norm(V2.(r2)) == 0.0

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

calculators = Any[]

# Basic lennard-jones calculator test
push!(calculators,
      (lennardjones(r0=rnn(:Al)),
       bulk(:Al, cubic=true, pbc=(true,false,false)) * (3,3,2) ) )

# PROBABLY NEED TO REVISIT THIS ONE
# # [3] JuLIP's EMT calculator
# at2 = set_pbc!( bulk(:Cu, cubic=true) * (2,2,2), (true,false,false) )
# set_positions!(at2, positions(at))
# emt2 = JuLIP.Potentials.EMTCalculator(at2)
# set_calculator!(at2, emt2)
# push!(calculators, (emt2, at2))

# Stillinger-Weber model
at3 = set_pbc!( bulk(:Si, cubic=true) * 2, (false, true, false) )
sw = StillingerWeber()
set_calculator!(at3, sw)
push!(calculators, (sw, at3))

# PairSitePotential (pair potential wrapped in a site potential)
at8 = set_pbc!( bulk(:Al, cubic=true), false ) * 2
pp = lennardjones(r0=rnn(:Al))
psp = SitePotential(pp)
if notCI
   push!(calculators, (psp, at8))
end

println("--------------------------------------------------")
println(" PairSitePotential Consistency test: ")
println("--------------------------------------------------")
println(" E_pp - E_psp = ", energy(pp, at8) - energy(psp, at8))
println(" |Frc_pp - Frc_psp| = ", maxnorm(forces(pp, at8) - forces(psp, at8)))
println("--------------------------------------------------")
@test abs(energy(pp, at8) - energy(psp, at8)) < 1e-11

# EAM Potential
at9 = set_pbc!( bulk(:Fe, cubic = true), false ) * (2,1,1)
eam = eam_Fe
push!(calculators, (eam, at9))

if eam_W4 != nothing
   # Another EAM Potential
   at10 = set_pbc!( bulk(:W, cubic = true), false ) * (2,1,2)
   eam4 = eam_W4
   push!(calculators, (eam4, at10))
end

# ========== Run the finite-difference tests for all calculators ============

println("============================================")
println("  Testing calculator implementations ")
println("============================================")
for (calc, at) in calculators
   println("--------------------------------")
   println(typeof(calc))
   @show length(at)
   println("--------------------------------")
   @test fdtest(calc, at, verbose=true)
end

# ========== Test correct implementation of site_energy ============
#            and of partial_energy

println("--------------------------------------------------")
println("Testing `site_energy` and `partial_energy` ...")
println("--------------------------------------------------")
at = bulk(:Si, pbc=true, cubic=true) * 3
sw = StillingerWeber()
atsm = bulk(:Si, pbc = true)
println("checking site energy identity . . .")
@test abs( JuLIP.Potentials.site_energy(sw, at, 1) - energy(sw, atsm) / 2 ) < 1e-10
rattle!(at, 0.01)
@test abs( energy(sw, at) - sum(site_energies(sw, at)) ) < 1e-10

println("fd test for site_energy")
# finite-difference test
set_constraint!(at, FixedCell(at))
f(x) = JuLIP.Potentials.site_energy(sw, set_dofs!(at, x), 1)
df(x) = (JuLIP.Potentials.site_energy_d(sw, set_dofs!(at, x), 1) |> mat)[:]
@test fdtest(f, df, dofs(at); verbose=true)

println("fd test for partial_energy")
Idom = [2,4,10]
f(x) = JuLIP.Potentials.partial_energy(sw, set_dofs!(at, x), Idom)
df(x) = (JuLIP.Potentials.partial_energy_d(sw, set_dofs!(at, x), Idom) |> mat)[:]
@test fdtest(f, df, dofs(at); verbose=true)

#  TODO: put these back in at some point . . .
# ================================================================
# # [5] a simple FDPotential
# @pot type FDPot <: FDPotential end
# fdpot(r) = exp(-0.3*r) * JuLIP.Potentials.cutsw(r, 4.0, 1.0)
# JuLIP.Potentials.ad_evaluate{T<:Real}(pot::FDPot, R::Matrix{T}) =
#                sum( fdpot(Base.LinAlg.vecnorm2(R[:,i])) for i = 1:size(R,2) )
# JuLIP.cutoff(::FDPot) = 4.0
# at5 = set_pbc!(bulk("Si") * (3,3,1), false)
# push!(calculators, (FDPot(), at5))
#
# # [6] a simple FDPotential
# @pot type FDPot_r <: FDPotential_r end
# JuLIP.Potentials.ad_evaluate{T<:Real}(pot::FDPot_r, r::Vector{T}) = sum( fdpot.(r) )
# JuLIP.cutoff(::FDPot_r) = 4.0
# at6 = set_pbc!(bulk("Si") * (3,3,1), false)
# push!(calculators, (FDPot_r(), at6))
#
# # [7] a simple RDPotential
# @pot type RDPot_r <: RDPotential_r end
# JuLIP.Potentials.ad_evaluate{T<:Real}(pot::RDPot_r, r::Vector{T}) = sum( fdpot.(r) )
# JuLIP.cutoff(::RDPot_r) = 4.0
# at7 = set_pbc!(bulk("Si") * (3,3,1), false)
# push!(calculators, (RDPot_r(), at7))
# ================================================================

ext: jl


text:
function push_to_abc_log!(abc_log::abc_log_type, plan::abc_pmc_plan_type, theta_star::Array{Float64,1}, ss_star, dist_star::Float64)
         if plan.save_params
           push!(abc_log.theta, theta_star)
         end
         if plan.save_summary_stats
           push!(abc_log.ss, ss_star)
         end
         if plan.save_distances
           push!(abc_log.dist, dist_star)
         end
end

function append_to_abc_log!(abc_log::abc_log_type, plan::abc_pmc_plan_type, theta_star::Array{Array{Float64,1},1}, ss_star::Array{Any,1}, dist_star::Array{Float64,1})
         if plan.save_params
           append!(abc_log.theta, theta_star)
         end
         if plan.save_summary_stats
           append!(abc_log.ss, ss_star)
         end
         if plan.save_distances
           append!(abc_log.dist, dist_star)
         end
end

function draw_indices(abc_log::ABC.abc_log_type, num_use::Integer, start_idx::Integer = 1, stop_idx::Integer = length(abc_log.dist))
  @assert num_use >=1
  num_available = 1+stop_idx-start_idx
  idx = collect(start_idx:min(start_idx+num_available-1,start_idx+num_use-1))
end

function draw_indicies_from_generations(accept_log::ABC.abc_log_type, reject_log::ABC.abc_log_type, generation_list::Vector{Int64}, num_use::Integer = 10)
 @assert length(generation_list) >= 1
 num_use_accepts = num_use
 num_use_rejects = num_use
 idx_accept = Int64[]
 idx_reject = Int64[]
 for generation in generation_list
    if (generation<1)||(generation>length(accept_log.generation_starts_at)) continue end
    @assert length(accept_log.generation_starts_at) == length(reject_log.generation_starts_at) >= generation >=1
    stop_idx_accept = generation+1 <= length(accept_log.generation_starts_at) ? accept_log.generation_starts_at[generation+1] : length(accept_log.dist)
    stop_idx_reject = generation+1 <= length(reject_log.generation_starts_at) ? reject_log.generation_starts_at[generation+1] : length(reject_log.dist)
    idx_accept_this_gen = draw_indices(accept_log,num_use_accepts,accept_log.generation_starts_at[generation],stop_idx_accept)
    num_use_rejects_adjusted = num_use_rejects + num_use_accepts - length(idx_accept_this_gen)
    idx_reject_this_gen = draw_indices(reject_log,num_use_rejects_adjusted,reject_log.generation_starts_at[generation],stop_idx_reject)
    append!(idx_accept,idx_accept_this_gen)
    append!(idx_reject,idx_reject_this_gen)
 end
 return (idx_accept,idx_reject)
end

function make_training_data(accept_log::ABC.abc_log_type, reject_log::ABC.abc_log_type, generation_list::Vector{Int64}, num_use::Integer = 10)
 @assert length(accept_log.theta) == length(accept_log.ss) == length(accept_log.dist) >= 1
 @assert length(reject_log.theta) == length(reject_log.ss) == length(reject_log.dist) >= 1
 @assert length(accept_log.theta[1]) >=1
 @assert length(reject_log.theta[1]) >=1
 @assert length(accept_log.ss[1]) == length(reject_log.ss[1])
 (idx_accept,idx_reject) = draw_indicies_from_generations(accept_log,reject_log,generation_list,num_use)

 num_use_total = length(idx_accept) + length(idx_reject)
 num_param = length(accept_log.theta[1])
 num_ss = length(accept_log.ss[1])
 x_train = Array{Float64}(num_param,num_use_total)
 y_train = Array{Float64}(num_ss,num_use_total)
 println("# idx_accept = ",idx_accept)
 println("# size(x_train)",size(x_train))
 println("# size(accept_log.theta[idx_accept]) = ", size( accept_log.theta[idx_accept]))
 for i in 1:length(idx_accept)
    x_train[:,i] = accept_log.theta[idx_accept[i]]
    y_train[:,i] = accept_log.ss[idx_accept[i]]
 end
 for i in 1:length(idx_reject)
    x_train[:,length(idx_accept)+i] = reject_log.theta[idx_reject[i]]
    y_train[:,length(idx_accept)+i] = reject_log.ss[idx_reject[i]]
 end
 return (x_train, y_train)
end

###=
function test_emu(pop_out::abc_population_type, theta::Array{Float64,1}, num_use::Integer = 100)
num_use = 100
(x_train,y_train) = ABC.make_training_data(pop_out.accept_log, pop_out.reject_log, [length(pop_out.accept_log.generation_starts_at)], num_use)
 sample_var_y_train = var(y_train,2)
 sigmasq_y_train = similar(y_train)
 for i in 1:size(y_train,2)
   sigmasq_y_train[:,i] = sample_var_y_train
 end
 emu = ABC.train_gp(x_train,y_train,sigmasq_y_train)
 function emulate_ss(theta::Array{Float64,2})
   predict_gp(emu,theta)
 end
 function emulate_ss(theta::Array{Float64,1})
   predict_gp(emu,reshape(theta,size(theta,1),1))
 end

 function emulator_output_to_ss_mean_stddev( emu_output::Tuple{Array{Float64,2},Array{Float64,3}} )
   m = emu_output[1]
   stddev = similar(emu_output[1])
   for i in 1:size(emu_output[1],1)
     stddev[i,:] = diag(emu_output[2][1,:,:])
   end
   return (m,stddev)
 end

  emu_out = emulate_ss(theta)
  ss_mean_stddev = emulator_output_to_ss_mean_stddev(emu_out)
  return (vec(ss_mean_stddev[1]),vec(ss_mean_stddev[2]))
end

function emu_accept_prob(ss_true::Array{Float64,1},emu_out::Tuple{Array{Float64,1},Array{Float64,1}}, epsilon::Real)
  mu = emu_out[1]
  sig = emu_out[2]
  prob_accept = 1.0
  for i in 1:length(ss_true)
    cdf_hi = cdf(Normal(mu[i],sig[i]),ss_true[i]+epsilon)
    cdf_lo = cdf(Normal(mu[i],sig[i]),ss_true[i]-epsilon)
    prob_accept *= cdf_hi-cdf_lo
  end
  return prob_accept
end

function emu_accept_reject_run_full_model(ss_true::Array{Float64,1},emu_out::Tuple{Array{Float64,1},Array{Float64,1}}, epsilon::Real; prob_accept_crit::Real = 0.999, prob_reject_crit::Real = 0.001)
  prob_accept = emu_accept_prob(ss_true,emu_out,epsilon)
  if prob_accept > prob_accept_crit
    return 1
  elseif  prob_accept < prob_reject_crit
    return -1
  else
    return 0
  end
end

#=
x = collect(linspace(0.0,0.1,200))
y=  map(x->emu_accept_reject(ss_true,emu_out,x),x)
plot(x,y)
=#


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

name = "x265"
version = v"3.5"

# Collection of sources required to build x265
sources = [
    # Note: x265 is only available as a Mercurial repository at
    # http://hg.videolan.org/x265/ They aren't publishing release tarballs
    # anymore.  This source is a git mirror on BitBucket.
    ArchiveSource("https://bitbucket.org/multicoreware/x265_git/get/Release_$(version.major).$(version.minor).tar.bz2",
                  "8692745fc050930e16dff725895c9201b3a9812467f570ada2f156f44b153c71"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/*x265*/
# Remove `-march` and `-mcpu` flags
SED_SCRIPTS=(-e  -e -- 's/-mcpu=native //g')
for CMAKE_FILE in source/CMakeLists.txt source/dynamicHDR10/CMakeLists.txt; do
    sed -i 's/add_definitions(-march=i686)//g' "${CMAKE_FILE}"
    sed -i 's/-mcpu=native //g' "${CMAKE_FILE}"
done
mkdir bld && cd bld
FLAGS=()
if [[ "${target}" == i686-* ]] || [[ "${target}" == aarch64-apple-darwin* ]]; then
    FLAGS+=(-DENABLE_ASSEMBLY=OFF)
fi
cmake -DCMAKE_INSTALL_PREFIX="${prefix}" \
    -DCMAKE_TOOLCHAIN_FILE="${CMAKE_TARGET_TOOLCHAIN}" \
    -DENABLE_PIC=ON \
    -DENABLE_SHARED=ON \
    ${FLAGS[@]} \
    ../source
make -j${nproc}
make install
# Remove the large static archive
rm  ${prefix}/lib/libx265.a
"""

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

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

# Dependencies that must be installed before this package can be built
dependencies = [
    # We need `nasm` for x86_64
    HostBuildDependency("NASM_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# This file is a part of RRRMC.jl. License is MIT: http://github.com/carlobaldassi/RRRMC.jl/LICENCE.md

module LEAliases

using ..LE
using ..Empty
using ..SK
using ..EA
using ..SAT
using ..PercLinear
using ..PercStep
using ..CommStep
using ..CommReLU
using ..CommQu

export Graph0LE, GraphSKLE, GraphEALE, GraphSATLE,
       GraphPercLinearLE, GraphPercStepLE, GraphCommStepLE,
       GraphCommReLULE, GraphCommQuLE

const Graph0LE{M,γT} = GraphLocalEntropy{M,γT,GraphEmpty}

# """
#     Graph0LE(...)
#
# TODO
# """
Graph0LE(Nk::Integer, M::Integer, γ::Float64, β::Float64) = GraphLocalEntropy(Nk, M, γ, β, GraphEmpty, Nk)

const GraphSKLE{M,γT} = GraphLocalEntropy{M,γT,GraphSK}

# """
#     GraphSKLE(N::Integer, M::Integer, γ::Float64, β::Float64) <: DoubleGraph
#
# TODO
# """
GraphSKLE(Nk::Integer, M::Integer, γ::Float64, β::Float64) = GraphLocalEntropy(Nk, M, γ, β, GraphSK, SK.gen_J(Nk))

const GraphEALE{M,γT,twoD} = GraphLocalEntropy{M,γT,GraphEANormal{twoD}}

# """
#     GraphEALE(L::Integer, D::Integer, M::Integer, γ::Float64, β::Float64) <: DoubleGraph
#
# TODO
# """
function GraphEALE(L::Integer, D::Integer, M::Integer, γ::Float64, β::Float64)
    D ≥ 1 || throw(ArgumentError("D must be ≥ 0, given: $D"))
    A = EA.gen_EA(L, D)
    N = length(A)
    J = EA.gen_J(Float64, N, A) do
        4*rand() - 2
    end

    GraphLocalEntropy(N, M, γ, β, GraphEANormal{2D}, L, A, J)
end

function GraphEALE(fname::AbstractString, M::Integer, γ::Float64, β::Float64)
    L, D, A, J = EA.gen_AJ(fname)
    N = length(A)
    @assert N == L^D
    GraphLocalEntropy(N, M, γ, β, GraphEANormal{2D}, L, A, J)
end

function GraphEALE(X::GraphEANormal{twoD}, M::Integer, γ::Float64, β::Float64) where {twoD}
    @assert iseven(twoD)
    D = twoD ÷ 2
    N = X.N
    L = round(Int, N^(1/D))
    @assert L^D == N
    GraphLocalEntropy(N, M, γ, β, GraphEANormal{twoD}, L, X.A, X.J)
end

const GraphSATLE{M,γT} = GraphLocalEntropy{M,γT,GraphSAT}

# """
#     GraphSATLE(N::Integer, K::Integer, α::Real, M::Integer, γ::Float64, β::Float64) <: DoubleGraph
#
# TODO
# """
function GraphSATLE(N::Integer, K::Integer, α::Real, M::Integer, γ::Float64, β::Float64)
    A, J = SAT.gen_randomKSAT(N, K, α)
    GraphLocalEntropy(N, M, γ, β, GraphSAT, N, A, J)
end

function GraphSATLE(X::GraphSAT, M::Integer, γ::Float64, β::Float64)
    N = X.N
    GraphLocalEntropy(N, M, γ, β, GraphSAT, N, X.A, X.J)
end

const GraphPercLinearLE{M,γT} = GraphLocalEntropy{M,γT,GraphPercLinear}

# """
#     GraphPercLinearLE(N::Integer, P::Integer, M::Integer, γ::Float64, β::Float64) <: DoubleGraph
#
# TODO
# """
function GraphPercLinearLE(N::Integer, P::Integer, M::Integer, γ::Float64, β::Float64)
    ξ, ξv = PercLinear.gen_ξ(N, P)
    GraphLocalEntropy(N, M, γ, β, GraphPercLinear, ξ, ξv)
end

function GraphPercLinearLE(X::GraphPercLinear, M::Integer, γ::Float64, β::Float64)
    GraphLocalEntropy(X.N, M, γ, β, GraphPercLinear, X.ξ, X.ξv)
end

const GraphPercStepLE{M,γT} = GraphLocalEntropy{M,γT,GraphPercStep}

# """
#     GraphPercStepLE(N::Integer, P::Integer, M::Integer, γ::Float64, β::Float64) <: DoubleGraph
#
# TODO
# """
function GraphPercStepLE(N::Integer, P::Integer, M::Integer, γ::Float64, β::Float64)
    ξ, ξv = PercStep.gen_ξ(N, P)
    GraphLocalEntropy(N, M, γ, β, GraphPercStep, ξ, ξv)
end

function GraphPercStepLE(X::GraphPercStep, M::Integer, γ::Float64, β::Float64)
    GraphLocalEntropy(X.N, M, γ, β, GraphPercStep, X.ξ, X.ξv)
end

const GraphCommStepLE{M,γT} = GraphLocalEntropy{M,γT,GraphCommStep}

# """
#     GraphCommStepLE(K1::Integer, K2::Integer, P::Integer, M::Integer, γ::Float64, β::Float64) <: DoubleGraph
#
# TODO
# """
function GraphCommStepLE(K1::Integer, K2::Integer, P::Integer, M::Integer, γ::Float64, β::Float64; fc::Bool = false)
    N = K1 * K2
    Kin = fc ? K1 : N
    ξ, ξv = CommStep.gen_ξ(Kin, P)
    if fc
        ξ = repeat(ξ, outer=(1,K2))
        ξv = [repeat(ξ1, K2) for ξ1 in ξv]
    end
    GraphLocalEntropy(N, M, γ, β, GraphCommStep, K2, ξ, ξv)
end

function GraphCommStepLE(X::GraphCommStep, M::Integer, γ::Float64, β::Float64)
    GraphLocalEntropy(X.N, M, γ, β, GraphCommStep, X.K2, X.ξ, X.ξv)
end

# """
#     GraphCommReLULE(K1::Integer, K2::Integer, P::Integer, M::Integer, γ::Float64, β::Float64) <: DoubleGraph
#
# TODO
# """
function GraphCommReLULE(K1::Integer, K2::Integer, P::Integer, M::Integer, γ::Float64, β::Float64; fc::Bool = false)
    N = K1 * K2
    Kin = fc ? K1 : N
    ξ, ξv, y = CommReLU.gen_ξ(Kin, P)
    if fc
        ξ = repeat(ξ, outer=(1,K2))
        ξv = [repeat(ξ1, K2) for ξ1 in ξv]
    end
    GraphLocalEntropy(N, M, γ, β, GraphCommReLU, K2, ξ, ξv, y)
end

function GraphCommReLULE(X::GraphCommReLU, M::Integer, γ::Float64, β::Float64)
    GraphLocalEntropy(X.N, M, γ, β, GraphCommReLU, X.K2, X.ξ, X.ξv, X.y)
end

# """
#     GraphCommQuLE(K1::Integer, K2::Integer, P::Integer, M::Integer, γ::Float64, β::Float64) <: DoubleGraph
#
# TODO
# """
function GraphCommQuLE(K1::Integer, K2::Integer, P::Integer, M::Integer, γ::Float64, β::Float64; fc::Bool = false)
    N = K1 * K2
    Kin = fc ? K1 : N
    ξ, ξv, y = CommReLU.gen_ξ(Kin, P)
    if fc
        ξ = repeat(ξ, outer=(1,K2))
        ξv = [repeat(ξ1, K2) for ξ1 in ξv]
    end
    GraphLocalEntropy(N, M, γ, β, GraphCommQu, K2, ξ, ξv, y)
end

function GraphCommQuLE(X::GraphCommQu, M::Integer, γ::Float64, β::Float64)
    GraphLocalEntropy(X.N, M, γ, β, GraphCommQu, X.K2, X.ξ, X.ξv, X.y)
end

end # module


text:
function is_func_call(x)
    if typof(x) === Call
        return true
    elseif typof(x) === WhereOpCall
        return is_func_call(x.args[1])
    elseif typof(x) === InvisBrackets
        return is_func_call(x.args[2])
    elseif typof(x) === UnaryOpCall
        return !(isoperator(x.args[1]) && (kindof(x.args[1]) === Tokens.EX_OR || kindof(x.args[1]) === Tokens.DECLARATION))
    elseif typof(x) === BinaryOpCall
        if issyntaxcall(x.args[2])
            if is_decl(x.args[2])
                return is_func_call(x.args[1])
            else
                return false
            end
        else
            true
        end
    else
        return false
    end
end

is_assignment(x) = typof(x) === BinaryOpCall && kindof(x.args[2]) === Tokens.EQ

# OPERATOR
is_exor(x) = isoperator(x) && kindof(x) == Tokens.EX_OR && x.dot == false
is_decl(x) = isoperator(x) && kindof(x) == Tokens.DECLARATION
is_issubt(x) = isoperator(x) && kindof(x) == Tokens.ISSUBTYPE
is_issupt(x) = isoperator(x) && kindof(x) == Tokens.ISSUPERTYPE
is_and(x) = isoperator(x) && kindof(x) == Tokens.AND && x.dot == false
is_not(x) = isoperator(x) && kindof(x) == Tokens.NOT && x.dot == false
is_plus(x) = isoperator(x) && kindof(x) == Tokens.PLUS && x.dot == false
is_minus(x) = isoperator(x) && kindof(x) == Tokens.MINUS && x.dot == false
is_star(x) = isoperator(x) && kindof(x) == Tokens.STAR && x.dot == false
is_eq(x) = isoperator(x) && kindof(x) == Tokens.EQ && x.dot == false
is_dot(x) = isoperator(x) && kindof(x) == Tokens.DOT
is_ddot(x) = isoperator(x) && kindof(x) == Tokens.DDOT
is_dddot(x) = isoperator(x) && kindof(x) == Tokens.DDDOT
is_pairarrow(x) = isoperator(x) && kindof(x) == Tokens.PAIR_ARROW && x.dot == false
is_in(x) = isoperator(x) && kindof(x) == Tokens.IN && x.dot == false
is_elof(x) = isoperator(x) && kindof(x) == Tokens.ELEMENT_OF && x.dot == false
is_colon(x) = isoperator(x) && kindof(x) == Tokens.COLON
is_prime(x) = isoperator(x) && kindof(x) == Tokens.PRIME
is_cond(x) = isoperator(x) && kindof(x) == Tokens.CONDITIONAL
is_where(x) = isoperator(x) && kindof(x) == Tokens.WHERE
is_anon_func(x) = isoperator(x) && kindof(x) == Tokens.ANON_FUNC

# PUNCTUATION
is_punc(x) = typof(x) === PUNCTUATION &&
    kindof(x) == Tokens.COMMA &&
    kindof(x) == Tokens.LPAREN &&
    kindof(x) == Tokens.RPAREN &&
    kindof(x) == Tokens.LBRACE &&
    kindof(x) == Tokens.RBRACE &&
    kindof(x) == Tokens.LSQUARE &&
    kindof(x) == Tokens.RSQUARE
is_comma(x) = ispunctuation(x) && kindof(x) == Tokens.COMMA
is_lparen(x) = ispunctuation(x) && kindof(x) == Tokens.LPAREN
is_rparen(x) = ispunctuation(x) && kindof(x) == Tokens.RPAREN
is_lbrace(x) = ispunctuation(x) && kindof(x) == Tokens.LBRACE
is_rbrace(x) = ispunctuation(x) && kindof(x) == Tokens.RBRACE
is_lsquare(x) = ispunctuation(x) && kindof(x) == Tokens.LSQUARE
is_rsquare(x) = ispunctuation(x) && kindof(x) == Tokens.RSQUARE

# KEYWORD
is_if(x) = iskw(x) && kindof(x) == Tokens.IF
is_module(x) = iskw(x) && kindof(x) == Tokens.MODULE
is_import(x) = iskw(x) && kindof(x) == Tokens.IMPORT
is_importall(x) = iskw(x) && kindof(x) == Tokens.IMPORTALL

# Literals
is_lit_string(x) = isliteral(x) && (kindof(x) == Tokens.STRING || kindof(x) == Tokens.TRIPLE_STRING)

function _arg_id(x)
    if typof(x) === IDENTIFIER
        return x
    elseif typof(x) === Quotenode
        return x.args[1]
    elseif typof(x) === Curly ||
        typof(x) === Kw ||
        typof(x) === WhereOpCall ||
           (typof(x) === UnaryOpCall && is_dddot(x.args[2])) ||
           (typof(x) === BinaryOpCall && is_decl(x.args[2])) ||
        return _arg_id(x.args[1])
    else
        return x
    end
end

function get_where_params(x, params = [])
    if typof(x) === WhereOpCall
        for i = 3:length(x.args)
            a = x.args[i]
            if !ispunctuation(a)
                param = rem_subtype(a)
                param = rem_curly(param)
                push!(params, str_value(param))
            end
        end
    end
    return params
end

function get_curly_params(x, params = [])
    if typof(x) == Curly
        for i = 2:length(x.args)
            a = x.args[i]
            if !ispunctuation(a)
                param = rem_subtype(a)
                push!(params, str_value(param))
            end
        end
    end
    return params
end

function get_sig_params(x, params = [])
    get_where_params(x, params)
    if typof(x) === WhereOpCall && typof(x.args[1]) === WhereOpCall
        get_where_params(x.args[1], params)
    end
    x = rem_where(x)
    x = rem_call(x)
    get_curly_params(x, params)
    return params
end

function rem_subtype(x)
    if typof(x) === BinaryOpCall && isoperator(x.args[2]) && kindof(x.args[2]) == Tokens.ISSUBTYPE
        return x.args[1]
    else
        return x
    end
end

function rem_decl(x)
    if typof(x) === BinaryOpCall && is_decl(x.args[2])
        return x.args[1]
    else
        return x
    end
end

function rem_curly(x)
    if typof(x) === Curly
        return x.args[1]
    else
        return x
    end
end

function rem_call(x)
    if typof(x) === Call
        return x.args[1]
    else
        return x
    end
end

function rem_where(x)
    if typof(x) === WhereOpCall
        return rem_where(x.args[1])
    else
        return x
    end
end

function rem_where_subtype(x)
    if typof(x) === WhereOpCall || typof(x) === BinaryOpCall && kindof(x.args[2]) === Tokens.ISSUBTYPE
        return rem_where_subtype(x.args[1])
    else
        return x
    end
end

function rem_where_decl(x)
    if typof(x) === WhereOpCall || typof(x) === BinaryOpCall && kindof(x.args[2]) === Tokens.DECLARATION
        return rem_where_decl(x.args[1])
    else
        return x
    end
end

function rem_invis(x)
    if typof(x) === InvisBrackets
        return x.args[2]
    else
        return x
    end
end

function rem_dddot(x)
    if typof(x) === UnaryOpCall && is_dddot(x.args[2])
        return x.args[1]
    else
        return x
    end
end

function rem_kw(x)
    if typof(x) === Kw
        return x.args[1]
    else
        return x
    end
end

# Definitions
function defines_function(x)
    if typof(x) === FunctionDef
        return true
    elseif typof(x) === BinaryOpCall
        if is_eq(x.args[2])
            sig = x.args[1]
            while true
                if typof(sig) === Call || typof(sig) === UnaryOpCall || (typof(sig) === BinaryOpCall && !(kindof(sig.args[2]) == Tokens.DOT || issyntaxcall(sig.args[2])))
                    return true
                elseif typof(sig) === BinaryOpCall && is_decl(sig.args[2]) || typof(sig) === WhereOpCall
                    sig = sig.args[1]
                else
                    return false
                end
            end
        end
        return false
    else
        return false
    end
end

defines_macro(x) = typof(x) == Macro
defines_datatype(x) = defines_struct(x) || defines_abstract(x) || defines_primitive(x)
defines_struct(x) = typof(x) === Struct || defines_mutable(x)
defines_mutable(x) = typof(x) === Mutable
defines_abstract(x) = typof(x) === Abstract
defines_primitive(x) = typof(x) === Primitive
defines_module(x) = typof(x) === ModuleH || typof(x) === BareModule
defines_anon_function(x) = typof(x) === BinaryOpCall && is_anon_func(x.args[2])

function has_sig(x)
    defines_datatype(x) || defines_function(x) || defines_macro(x) || defines_anon_function(x)
end

"""
    get_sig(x)

Returns the full signature of function, macro and datatype definitions.
Should only be called when has_sig(x) == true.
"""
function get_sig(x)
    if typof(x) === BinaryOpCall
        return x.args[1]
    elseif typof(x) === Struct ||
        typof(x) === FunctionDef ||
        typof(x) === Macro
        return x.args[2]
    elseif typof(x) === Mutable ||
           typof(x) === Abstract ||
           typof(x) === Primitive
        length(x.args) < 3 && error(x)
        return x.args[3]
    end
end

function get_name(x)
    if typof(x) === Struct || typof(x) === Mutable || typof(x) === Abstract || typof(x) === Primitive
        sig = get_sig(x)
        sig = rem_subtype(sig)
        sig = rem_where(sig)
        sig = rem_subtype(sig)
        sig = rem_curly(sig)
    elseif typof(x) === ModuleH || typof(x) === BareModule
        sig = x.args[2]
    elseif typof(x) === FunctionDef || typof(x) === Macro
        sig = get_sig(x)
        sig = rem_where(sig)
        sig = rem_decl(sig)
        sig = rem_call(sig)
        sig = rem_curly(sig)
        sig = rem_invis(sig)
        if typof(sig) === BinaryOpCall && kindof(sig.args[2]) == Tokens.DOT
            if length(sig.args) > 2 && sig.args[3].args isa Vector{EXPR} && length(sig.args[3].args) > 0
                sig = sig.args[3].args[1]
            end
        end
        return sig
        # return get_name(sig)
    elseif typof(x) === BinaryOpCall
        x.args === nothing && return x
        length(x.args) < 2 && return x
        if kindof(x.args[2]) == Tokens.DOT
            if length(x.args) > 2 && typof(x.args[3]) === Quotenode && x.args[3].args isa Vector{EXPR} && length(x.args[3].args) > 0
                return get_name(x.args[3].args[1])
            else
                return x
            end
        end
        sig = x.args[1]
        if typof(sig) === UnaryOpCall
            return get_name(sig.args[1])
        end
        sig = rem_where(sig)
        sig = rem_decl(sig)
        sig = rem_call(sig)
        sig = rem_curly(sig)
        sig = rem_invis(sig)
        return get_name(sig)
    else
        sig = x
        if typof(sig) === UnaryOpCall
            sig = sig.args[1]
        end
        sig = rem_where(sig)
        sig = rem_decl(sig)
        sig = rem_call(sig)
        sig = rem_curly(sig)
        sig = rem_invis(sig)
    end
end

function get_args(x)
    if typof(x) === IDENTIFIER
        return []
    elseif defines_anon_function(x) && !(typof(x.args[1]) === TupleH)
        arg = x.args[1]
        arg = rem_invis(arg)
        arg = get_arg_name(arg)
        return [arg]
    elseif typof(x) === TupleH
        args = []
        for i = 2:length(x.args)
            arg = x.args[i]
            ispunctuation(arg) && continue
            typof(arg) === Parameters && continue
            arg_name = get_arg_name(arg)
            push!(args, arg_name)
        end
        return args
    elseif typof(x) === Do
        args = []
        for i = 1:length(x.args[3].args)
            arg = x.args[3].args[i]
            ispunctuation(arg) && continue
            typof(arg) === Parameters && continue
            arg_name = get_arg_name(arg)
            push!(args, arg_name)
        end
        return args
    elseif typof(x) === Call
        args = []
        sig = rem_where(x)
        sig = rem_decl(sig)
        if typof(sig) === Call
            for i = 2:length(sig.args)
                arg = sig.args[i]
                ispunctuation(arg) && continue
                if typof(arg) === Parameters
                    for j = 1:length(arg.args)
                        parg = arg.args[j]
                        ispunctuation(parg) && continue
                        parg_name = get_arg_name(parg)
                        push!(args, parg_name)
                    end
                else
                    arg_name = get_arg_name(arg)
                    push!(args, arg_name)
                end
            end
        else
            error("not sig: $sig")
        end
        return args
    elseif typof(x) === Struct
        args = []
        for arg in x.args[3]
            if !defines_function(arg)
                arg = rem_decl(arg)
                push!(args, arg)
            end
        end
        return args
    elseif typof(x) === Mutable
        args = []
        for arg in x.args[4]
            if !defines_function(arg)
                arg = rem_decl(arg)
                push!(args, arg)
            end
        end
        return args
    elseif typof(x) === Flatten
        return get_args(x.args[1])
    elseif typof(x) === Generator || typof(x) === Flatten
        args = []
        if typof(x.args[1]) === Flatten || typof(x.args[1]) === Generator
            append!(args, get_args(x.args[1]))
        end

        if typof(x.args[3]) === Filter
            return get_args(x.args[3])
        else
            for i = 3:length(x.args)
                arg = x.args[i]
                if is_range(arg)
                    arg = rem_decl(arg.args[1])
                    arg = flatten_tuple(arg)
                    arg = rem_decl.(arg)
                    append!(args, arg)
                end
            end
            return args
        end
    else
        sig = get_sig(x)
        sig = rem_where(sig)
        sig = rem_decl(sig)
        return get_args(sig)
    end
end

function get_arg_name(arg)
    arg = rem_kw(arg)
    arg = rem_dddot(arg)
    arg = rem_where(arg)
    arg = rem_decl(arg)
    arg = rem_subtype(arg)
    arg = rem_curly(arg)
    arg = rem_invis(arg)
end

function get_arg_type(arg)
    if typof(arg) === BinaryOpCall && is_decl(arg.args[2])
        return Expr(arg.args[3])
    else
        return :Any
    end
end

get_body(x) = typof(x) === Mutable ? x.args[4] : x.args[3]

function flatten_tuple(x, out = [])
    if typof(x) === TupleH
        for arg in x
            ispunctuation(arg) && continue
            flatten_tuple(arg, out)
        end
    elseif typof(x) === InvisBrackets
        return flatten_tuple(x.args[2], out)
    else
        push!(out, x)
    end
    return out
end

"""
    get_id(x)

Get the IDENTIFIER name of a variable, possibly in the presence of
type declaration operators.
"""
function get_id(x)
    if typof(x) === BinaryOpCall && (is_issubt(x.args[2]) || is_decl(x.args[2])) ||
        (typof(x) === UnaryOpCall && is_dddot(x.args[2])) ||
        typof(x) === WhereOpCall ||
        typof(x) === Curly
        return get_id(x.args[1])
    elseif typof(x) === InvisBrackets
        return get_id(x.args[2])
    else
        return x
    end
end

# """
#     get_t(x)

# Basic inference in the presence of type declarations.
# """
# get_t(x) = :Any
# function get_t(x::BinaryOpCall)
#     if is_decl(x.args[2])
#         return Expr(x.args[3])
#     else
#         return :Any
#     end
# end

# infer_t(x) = :Any
# function infer_t(x::LITERAL)
#     if kindof(x) == Tokens.INTEGER
#         return :Int
#     elseif kindof(x) == Tokens.FLOAT
#         return :Float64
#     elseif kindof(x) == Tokens.STRING
#         return :String
#     elseif kindof(x) == Tokens.TRIPLE_STRING
#         return :String
#     elseif kindof(x) == Tokens.CHAR
#         return :Char
#     elseif kindof(x) == Tokens.TRUE || kindof(x) == Tokens.FALSE
#         return :Bool
#     elseif kindof(x) == Tokens.CMD
#         return :Cmd
#     end
# end

# infer_t(x::EXPR{Vect}) = :(Array{Any,1})
# infer_t(x::EXPR{Vcat}) = :(Array{Any,N})
# infer_t(x::EXPR{TypedVcat}) = :(Array{$(Expr(x.args[1])),N})
# infer_t(x::EXPR{Hcat}) = :(Array{Any,2})
# infer_t(x::EXPR{TypedHcat}) = :(Array{$(Expr(x.args[1])),2})
# infer_t(x::EXPR{Quote}) = :Expr
# infer_t(x::EXPR{StringH}) = :String
# infer_t(x::EXPR{Quotenode}) = :QuoteNode

"""
    contributes_scope(x)
Checks whether the body of `x` is included in the toplevel namespace.
"""
contributes_scope(x) = typof(x) in (FileH, Begin, Block, Const, Global, Local, If, MacroCall, TopLevel)


text: struct S a; a end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
module HepApi

import HTTP: URI, HTTP.get
using HTTP

using JSON

export search_inspirehep

const QueryParameterDict = Dict{AbstractString, Union{Integer, AbstractString}}

struct Search
	record_type::AbstractString
	query::QueryParameterDict
end

function URI(sch::Search)
	path = "/api/" * sch.record_type
	uri = HTTP.URI(;
								 scheme="https", host="inspirehep.net",
								 path=path, 	 query=sch.query
								 )
	return uri
end

get(sch::Search) = HTTP.get(URI(sch))

function check_if_successful_resp(resp::HTTP.Messages.Response)
	if resp.status >= 400
		e = ErrorException("Failed to get HEP data. The response code is $(resp.status)")
		throw(e)
	end
end

function resptojson(resp::HTTP.Messages.Response)
	return JSON.parse(String(resp.body))
end

"""
`HepApi.search_inspirehep(record_type::AbstractString, qpd::QueryParameterDict)`

Uses the inspirehep REST api to download a JSON from inspirehep.net.
The [Inspirehep Github Page](https://github.com/inspirehep/rest-api-doc) is the reference of this code.

e.g. (json data of the top 5 most cited papers on inpsirehep.net)

```
lit_json = search_inspirehep("literature", Dict("sort"=>"mostcited", "size"=>5))
println(lit)
```

- **record_type** are base categories for a search and can take the following values: **literature**, **authors**, **conferences**, **seminars**, **journals**, **jobs**, **experiments**, and **data**

- **q** is the text to search by inpsirehep's search syntax can be used.

- **sort** is the option by which the JSON is ordered. (allowed sort field TBD)

- **page** the page number returned.

- **fields** a fiter for certain metadata fields.

"""
function search_inspirehep(record_type::AbstractString, qpd)
	sch = Search(record_type, qpd)
	resp = get(sch)
	check_if_successful_resp(resp)
	return resptojson(resp)
end

end # module


text:
# First millisecond of 2015.
const DISCORD_EPOCH = 1420070400000

# Discord's form of ID.
const Snowflake = UInt64

snowflake(s::Integer) = Snowflake(s)
snowflake(s::AbstractString) = parse(Snowflake, s)

# TODO: Put these in helpers?
snowflake2datetime(s::Snowflake) = unix2datetime(((s >> 22) + DISCORD_EPOCH) / 1000)
worker_id(s::Snowflake) = (s & 0x3e0000) >> 17
process_id(s::Snowflake) = (s & 0x1f000) >> 12
increment(s::Snowflake) = s & 0xfff

# Discord sends both Unix and ISO timestamps.
datetime(s::Int64) = unix2datetime(s / 1000)
datetime(s::AbstractString) = DateTime(replace(s, "+" => ".000+")[1:23], ISODateTimeFormat)
datetime(d::DateTime) = d

# Lower to something JSON-compatible.
lowered(x::DateTime) = round(Int, datetime2unix(x))
lowered(x::Union{Integer, Bool}) = x
lowered(x::Vector) = lowered.(x)
lowered(x::Nothing) = nothing
lowered(x) = JSON.lower(x)

# Define JSON.lower for a type.
macro lower(T)
    if supertype(eval(T)) <: Enum{<:Integer}
        quote
            JSON.lower(x::$T) = Int(x)
        end
    else
        quote
            function JSON.lower(x::$T)
                d = Dict{Symbol, Any}()

                for f in fieldnames($T)
                    v = getfield(x, f)
                    if !ismissing(v)
                        d[f] = lowered(v)
                    end
                end

                return d
            end
        end
    end
end

# Define Base.merge for a type.
macro merge(T)
    quote
        function Base.merge(a::$T, b::$T)
            vals = map(fieldnames($T)) do f
                va = getfield(a, f)
                vb = getfield(b, f)
                ismissing(vb) ? va : vb
            end

            return $T(vals...)
        end
        Base.merge(::Missing, x::$T) = x
        Base.merge(x::$T, ::Missing) = x
    end
end

# Compute the expression needed to extract field k from keywords.
field(k::QuoteNode, ::Type{Any}) = :(kwargs[$k])
field(k::QuoteNode, ::Type{Snowflake}) = :(snowflake(kwargs[$k]))
field(k::QuoteNode, ::Type{DateTime}) = :(datetime(kwargs[$k]))
field(k::QuoteNode, ::Type{T}) where T = :($T(kwargs[$k]))
field(k::QuoteNode, ::Type{Vector{Snowflake}}) = :(snowflake.(kwargs[$k]))
field(k::QuoteNode, ::Type{Vector{DateTime}}) = :(datetime.(kwargs[$k]))
field(k::QuoteNode, ::Type{Vector{T}}) where T = :($T.(kwargs[$k]))
function field(k::QuoteNode, ::Type{T}) where T <: Enum
    return :(kwargs[$k] isa Integer ? $T(Int(kwargs[$k])) :
             kwargs[$k] isa $T ? kwargs[$k] : $T(kwargs[$k]))
end
function field(k::QuoteNode, ::Type{Optional{T}}) where T
    return :(haskey(kwargs, $k) ? $(field(k, T)) : missing)
end
function field(k::QuoteNode, ::Type{Nullable{T}}) where T
    return :(kwargs[$k] === nothing ? nothing : $(field(k, T)))
end
function field(k::QuoteNode, ::Type{OptionalNullable{T}}) where T
    return :(haskey(kwargs, $k) ? $(field(k, Nullable{T})) : missing)
end

# Define constructors from keyword arguments and a Dict for a type.
macro constructors(T)
    TT = eval(T)
    args = map(f -> field(QuoteNode(f), fieldtype(TT, f)), fieldnames(TT))

    quote
        $(esc(T))(; kwargs...) = $(esc(T))($(args...))
        $(esc(T))(d::Dict{Symbol, Any}) = $(esc(T))(; d...)
        $(esc(T))(x::$(esc(T))) = x
    end
end

# Export all instances of an enum.
macro exportenum(T)
    TT = eval(T)
    quote
        $(map(x -> :(export $(Symbol(x))), instances(TT))...)
    end
end

# Format a type for a docstring.
function doctype(s::String)
    s = replace(s, "UInt64" => "Snowflake")
    s = replace(s, string(Int) => "Int")
    s = replace(s, "Discord." => "")
    s = replace(s, "Dates." => "")
    m = match(r"Array{([^{}]+),1}", s)
    m === nothing || (s = replace(s, m.match => "Vector{$(m.captures[1])}"))
    m = match(r"Union{Missing, Nothing, (.+)}", s)
    m === nothing || return replace(s, m.match => "OptionalNullable{$(m.captures[1])}")
    m = match(r"Union{Missing, (.+)}", s)
    m === nothing || return replace(s, m.match => "Optional{$(m.captures[1])}")
    m = match(r"Union{Nothing, (.+)}", s)
    m === nothing || return replace(s, m.match => "Nullable{$(m.captures[1])}")
    return s
end

# Update a type's docstring with field names and types.
macro fielddoc(T)
    TT = eval(T)
    fields = filter(n -> !startswith(string(n), "djl_"), collect(fieldnames(TT)))
    ns = collect(string.(fields))
    width = maximum(length, ns)
    map!(n -> rpad(n, width), ns, ns)
    ts = collect(map(f -> string(fieldtype(TT, f)), fields))
    map!(doctype, ts, ts)
    docs = join(map(t -> "$(t[1]) :: $(t[2])", zip(ns, ts)), "\n")

    quote
        doc = string(@doc $T)
        docstring = doc * "\n## Fields\n\n```\n" * $docs * "\n```\n"

        Base.CoreLogging.with_logger(Base.CoreLogging.NullLogger()) do
            @doc docstring $T
        end
    end
end

# Produce a random string.
randstring() = String(filter(!ispunct, map(i -> Char(rand(48:122)), 1:rand(1:20))))

# Produce a randomized value of a type.
mock(::Type{Bool}; kwargs...) = rand(Bool)
mock(::Type{DateTime}; kwargs...) = now()
mock(::Type{AbstractString}; kwargs...) = randstring()
mock(::Type{Dict{Symbol, Any}}; kwargs...) = Dict(:a => mock(String), :b => mock(Int))
mock(::Type{T}; kwargs...) where T <: AbstractString = T(randstring())
mock(::Type{T}; kwargs...) where T <: Integer = abs(rand(T))
mock(::Type{T}; kwargs...) where T <: Enum = instances(T)[rand(1:length(instances(T)))]
mock(::Type{Vector{T}}; kwargs...) where T = map(i -> mock(T; kwargs...), 1:rand(1:10))
mock(::Type{Set{T}}; kwargs...) where T = Set(map(i -> mock(T; kwargs...), 1:rand(1:10)))
mock(::Type{Optional{T}}; kwargs...) where T = mock(T; kwargs...)
mock(::Type{Nullable{T}}; kwargs...) where T = mock(T; kwargs...)
mock(::Type{OptionalNullable{T}}; kwargs...) where T = mock(T; kwargs...)

# Define a mock method for a type.
macro mock(T)
    quote
        function $(esc(:mock))(::Type{$T}; kwargs...)
            names = fieldnames($(esc(T)))
            types = map(TT -> fieldtype($(esc(T)), TT), names)
            args = Vector{Any}(undef, length(names))
            for (i, (n, t)) in enumerate(zip(names, types))
                args[i] = haskey(kwargs, n) ? kwargs[n] : mock(t; kwargs...)
            end
            return $(esc(T))(args...)
        end
    end
end

# Apply the above macros to a type.
macro boilerplate(T, exs...)
    macros = map(e -> e.value, exs)

    quote
        @static if :constructors in $macros
            @constructors $T
        end
        @static if :docs in $macros
            @fielddoc $T
        end
        @static if :export in $macros
            @exportenum $T
        end
        @static if :lower in $macros
            @lower $T
        end
        @static if :merge in $macros
            @merge $T
        end
        @static if :mock in $macros
            @mock $T
        end
    end
end

include("overwrite.jl")
include("role.jl")
include("guild_embed.jl")
include("attachment.jl")
include("voice_region.jl")
include("activity.jl")
include("embed.jl")
include("user.jl")
include("ban.jl")
include("integration.jl")
include("connection.jl")
include("emoji.jl")
include("reaction.jl")
include("presence.jl")
include("channel.jl")
include("webhook.jl")
include("invite_metadata.jl")
include("member.jl")
include("voice_state.jl")
include("message.jl")
include("guild.jl")
include("invite.jl")
include("audit_log.jl")


text:
using NeuroEvolution
using Cambrian
using JSON
include("xor.jl")

cfg = NeuroEvolution.get_config("../cfg/test.yaml")
id = "011fc676-b63a-438a-a19a-4c9fc91a052b"

@testset "Save / Load NEATIndividual" begin
    ind = NEATIndividual(cfg)
    build!(ind)
    s = string(ind)
    new_ind = NEATIndiv(cfg, s)

    @test typeof(new_ind) == NEATIndividual
    build!(new_ind)
end

@testset "Save / Load HyperNEATIndividual" begin
    hn_cfg = NeuroEvolution.get_config("../cfg/hyperneat.yaml")
    ind = HyperNEATIndividual(hn_cfg)
    build!(ind)
    s = string(ind)
    new_ind = NEATIndiv(hn_cfg, s)

    @test typeof(new_ind) == HyperNEATIndividual
    build!(new_ind)
end

@testset "Save population" begin
    e = NEAT(cfg, fitness_xor, cfg["n_in"])
    id = e.id
    step!(e)
    Cambrian.save_gen(e)
    l = readdir("gens/$(e.id)/0001")
    @test length(l) > 0
end

@testset "Load population" begin
    println("Load test deactivated for external testing purposes")
    # e = NEAT(cfg, fitness_xor, cfg["n_in"])
    # load_gen!(e, "$id/0001")
    # @test all(getfield.(e.population, :fitness) .!= [-Inf])
    # step!(e)
end


text:
# --------------------------------------------------------------------------------------------------
# Functions and types for general filters

struct Filter
	coefficients::OffsetArrays.OffsetVector{Float64, Vector{Float64}}
    support::AbstractRange{Int64}

    function Filter(coefficients::AbstractVector{Float64}, support::AbstractRange{Int64})
        new(OffsetArrays.OffsetVector(coefficients, support), support)
    end
end

coefficients(h::Filter) = h.coefficients
support(h::Filter) = h.support
Base.length(h::Filter) = h |> coefficients |> length

function support_boundaries(h::Filter)
    filter_support = support(h)

    return filter_support[1], filter_support[end]
end

function Base.show(io::IO, h::Filter)
	show(io, coefficients(h))
end

# TODO: How to handle indices with more than one integer outside support?
function Base.getindex(h::Filter, idx::Integer)
	if checkbounds(Bool, coefficients(h), idx)
		return coefficients(h)[idx]
	else
		return 0.0
	end
end

# --------------------------------------------------------------------------------------------------
# Functions and types for interacting with interior filters

struct InteriorFilter
	vanishing_moments::Int64
	filter::Filter

	function InteriorFilter(p, filter)
        if p < 0
			throw(DomainError(p, "Not a valid number of vanishing moments"))
        end

        new(p, filter)
	end
end

filter(h::InteriorFilter) = h.filter
vanishing_moments(h::InteriorFilter) = h.vanishing_moments
coefficients(h::InteriorFilter) = h |> filter |> coefficients
support(h::InteriorFilter) = h |> filter |> support
support_boundaries(h::InteriorFilter) = h |> filter |> support_boundaries
Base.length(h::InteriorFilter) = h |> filter |> length

Base.getindex(h::InteriorFilter, idx) = filter(h)[idx]

function Base.show(io::IO, h::InteriorFilter)
    left, right = support_boundaries(h)
	println(io, "Filter for Daubechies ", vanishing_moments(h),
			" scaling function on [", left, ", ", right, "]:")
    show(io, filter(h))
end

"""
	interior_filter(p::Int[, phase=:symmlet])

Filter for Daubechies scaling function with `p` vanishing moments and `phase` being `:symmlet` or
`:min` (minimum phase).
"""
function interior_filter(p::Integer, phase::Symbol=:symmlet)
	if p < 1
		throw(DomainError(p, "Interior filter must have at least 1 vanishing moment"))
	end

    return InteriorFilter(p, interior_filter(p, Val{phase}))
end

function interior_filter(p::Integer, ::Type{Val{:symmlet}})
    if !(2 <= p <= 8)
        throw(DomainError(p, "symmlet filter should have between 2 and 8 vanishing moments"))
    end

    Filter(INTERIOR_COEFFICIENTS[p], -p+1:p)
end

function interior_filter(p::Integer, ::Type{Val{:min}})
    coefficients = Wavelets.wavelet(Wavelets.WT.Daubechies{p}()).qmf
    Filter(coefficients, 0:2*p-1)
end

# --------------------------------------------------------------------------------------------------
# Functions and types for interacting with boundary filters

struct BoundaryFilters
    side::Sides
    vanishing_moments::Int64
    filters::Vector{Filter}

    function BoundaryFilters(side, p, filters)
        if !(2 <= p <= 8)
            throw(DomainError(p, "Vanishing moments should be between 2 and 8"))
        end

        new(side, p, filters)
    end
end

side(B::BoundaryFilters) = B.side
vanishing_moments(B::BoundaryFilters) = B.vanishing_moments
filters(B::BoundaryFilters) = B.filters

function Base.getindex(B::BoundaryFilters, idx::Integer)
    filters(B)[idx + 1]
end

function boundary_filters(p::Integer, side::Sides)
    if !(2 <= p <= 8)
        throw(DomainError(p, "Vanishing moments should be between 2 and 8"))
    end

    supports = [0:(p + 2*k) for k in 0:(p - 1)]

    if side == LEFT
        coefficients = LEFT_SCALING_COEFFICIENTS[p]
    else
        coefficients = RIGHT_SCALING_COEFFICIENTS[p]
    end

    filters = map(Filter, coefficients, supports)

    return BoundaryFilters(side, p, filters)
end


text:
using Printf
using LinearAlgebra
using DelimitedFiles
using Quaternions

#EXERCISE MADE BY ERIK MARTÍN GARZÓN & DAVID BOCES OBIS
#-------------------------------------------------------------------- EXERCISE 2.2

V3 = [0; 0; 1] #vector of R3 we want to project in R2

function to_2d(V3)
    #we calculate de matrix (see the first exercise)

    a = 0.5*cosd(42)
    b = 0.5*sind(42)
    c = cosd(7)
    d = sind(7)
    e = 0
    f = 1

    M = [-a c e;
        -b -d f]
     #matrix found in the first exercise

    V2 = M*V3 #operation to transform V3 of R3 in V2 in R2

    # println("The vector ", V3, " representated in R2 is: ", V2)
    # println()

    return V2
end

#-------------------------------------------------------------------- EXERCISE 2.3
function rotate_phi_z(phi, z, V)

    global z_norm
    global z_module

    z_module = sqrt(z[1]*z[1]+z[2]*z[2]+z[3]*z[3])

    z_norm = z/z_module

    Z = [0.0 -z_norm[3] z_norm[2];
         z_norm[3] 0.0 -z_norm[1];
        -z_norm[2] z_norm[1] 0.0]

    #formula de Rodriges to find "R", the rotation matrix
    R = I + sind(phi)*Z + (1-cosd(phi))*(Z^2)
    #just above, we multiply pi/180 to the angle (in degrees) to obtain the angle in radians
    #so julia can do the operation we want

    #next, we will check if R is a rotation matrix: R*R'=I and det(R)=1
    #For this comprobation, we needed to round R*R' and det(R). That's because much times
    #the result was 0.999999 and not 1
    if(phi==0)
        return V
    end
    if (z != [0.0;0.0;0.0])
        if (round.(R*R') == I && round(det(R)) == 1)
            #we find the rotated matrix:
            U = R*V
            # print("V rotated is = ", U)
            # println()
            return U
        else

            print("R = ", R, "\nis NOT a rotation matrix\n")
            println("R*R'= ", round.(R*R'))
            println("det(R)= ", round(det(R)))

        end
    else

            return [0.0 0.0 0.0; 0.0 0.0 0.0;0.0 0.0 0.0]
    end
end

#-------------------------------------------------------------------- EXERCISE 3.1
global axis_norm

function axis_angle_to_mat(axis, angle)
    #first, we normalize the axis
    modul = sqrt(axis[1]*axis[1] + axis[2]*axis[2] + axis[3]*axis[3])

    global axis_norm

    round(modul)

    if (modul == 0)
        return [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
    else
        axis_norm = axis/modul

        Z = [0 -axis_norm[3] axis_norm[2];
             axis_norm[3] 0 -axis_norm[1];
            -axis_norm[2] axis_norm[1] 0]

        #formula Rodrigues
        R = I + sind(angle)*Z + (1-cosd(angle))*(Z^2)
        #just above, we multiply pi/180 to the angle (in degrees) to obtain the angle in radians
        #so julia can do the operation we want

        #next, we will check if R is a rotation matrix: R*R'=I and det(R)=1
        #For this comprobation, we needed to round R*R' and det(R). That's because much times
        #the result was 0.999999 and not 1
        if (round.(R*R')== I && round(det(R))==1)
            #we return the rotation matrix
            return R
        else
            print("R = ", R, "\nis NOT a rotation matrix\n")
            println("R*R'= ", round.(R*R'))
            println("det(R)= ", round(det(R)))
        end
    end
end

function mat_to_axis_angle(mat)
    #first, we find the angle:
    angle = acosd((tr(mat)-1)/2)

    #then, we find the axis:
    V = (mat-mat')/(2*sind(angle))
    v = [V[3,2]; V[1,3]; V[2,1]]

    #we return the angle in degrees
    return angle, v
end

function axis_angle_to_quat(axis, angle)
    global q2
    #first, we normalize the axis
    modul = sqrt(axis[1]*axis[1] + axis[2]*axis[2] + axis[3]*axis[3])

# if (modul==0)
#         print("\n QUATERNION MODULE IS 0")
#         return 0
#     else
        global axis_norm = axis/modul

        q2 = quat(cosd(angle/2), sind(angle/2)*axis_norm[1], sind(angle/2)*axis_norm[2], sind(angle/2)*axis_norm[3])

        return q2
    # end
end

function quat_to_axis_angle(q3)

    norm = sqrt(q3.s*q3.s + q3.v1*q3.v1 + q3.v2*q3.v2 + q3.v3*q3.v3)

    global axis_norm = [0.0;0.0;0.0]

    if(norm == 0)
        angle3 = 0

        axis_norm[1] = 0
        axis_norm[2] = 0
        axis_norm[3] = 0

    else
        angle3 = 2*acosd(q3.s/norm)
    end
    if (angle3==0)
        # println("Angle of the quaternion = 0  Returning 0 axis")
        return angle3, axis_norm
    else
        axis_norm[1] = (q3.v1/norm)/sind(angle3/2)
        axis_norm[2] = (q3.v2/norm)/sind(angle3/2)
        axis_norm[3] = (q3.v3/norm)/sind(angle3/2)
    end
    return angle3, axis_norm
end

function quat_to_mat(q)
    global matQuat = [((q.s)^2 + (q.v1)^2 - (q.v2)^2 - (q.v3)^2) (2*(q.v1)*(q.v2) - 2*(q.s)*(q.v3)) (2*(q.v1)*(q.v3) + 2*(q.s)*(q.v2));
            (2*(q.v1)*(q.v2) + 2*(q.s)*(q.v3)) ((q.s)^2 - (q.v1)^2 + (q.v2)^2 - (q.v3)^2) (2*(q.v2)*(q.v3) - 2*(q.s)*(q.v1));
            (2*(q.v1)*(q.v3) - 2*(q.s)*(q.v2)) (2*(q.v2)*(q.v3) + 2*(q.s)*(q.v1)) ((q.s)^2 - (q.v1)^2 - (q.v2)^2 + (q.v3)^2)]

    return matQuat
end

function mat_to_quat(mat)
    #first, we find the angle:
    angle2 = acosd((tr(mat)-1)/2)

    if (angle2 != 0)
        #then, we find the axis:
        V = (mat-mat')/(2*sind(angle2))
        v = [V[3,2]; V[1,3]; V[2,1]]

        q3 = quat(cosd(angle2/2), sind(angle2/2)*v[1], sind(angle2/2)*v[2], sind(angle2/2)*v[3])

        return q3
    else
        return "THIS ROTATION MATRIX [IDENTITY MATRIX] HAS NO AXIS. ITS ROTATION ANGLE IS 0 + 360*K [BEING K AN ENTER]. FOR THE RESULT THIS MATRIX HAS NO QUATERNION FORM."
    end
end

#-------------------------------------------------------------------- EXERCISE 4.3
function scale_and_translation(v3d, scale, point)

# println("V3D ", v3d)

    v3d = v3d*scale

    # println("V3D SCALED ", v3d)

    #movement = [0.0;0.0]
    v3dResult = [0.0;0.0;0.0]

    # movement[1] = point[1] - v2d[1]
    # movement[2] = point[2] - v2d[2]

    v3dResult[1] = v3d[1] + point[1]
    v3dResult[2] = v3d[2] + point[2]
    v3dResult[3] = v3d[3] + point[3]

    return v3dResult

end


text:
using Distributed; addprocs()
@everywhere begin
    using Pkg
    Pkg.activate(Base.source_dir())
end

using DiscreteChoiceModels, Dagger, Tables, CSV, Infiltrator

function main()
    data = DTable(Tables.rowtable(CSV.File("test/data/biogeme_swissmetro.dat", delim='\t')), 2000)

    # filter to wanted cases
    data = filter(data) do row
        ((row.PURPOSE == 1) | (row.PURPOSE == 3)) && row.CHOICE != 0
    end

    data = map(data) do r
        merge((
            avtr=(r.TRAIN_AV == 1) && (r.SP != 0),
            avsm=r.SM_AV == 1,
            avcar=(r.CAR_AV == 1) && (r.SP != 0)
        ), pairs(r))
    end

    @time multinomial_logit(
        @utility(begin
            1 ~ αtrain + βtravel_time * TRAIN_TT / 100 + βcost * (TRAIN_CO * (GA == 0)) / 100
            2 ~ αswissmetro + βtravel_time * SM_TT / 100 + βcost * SM_CO * (GA == 0) / 100
            3 ~ αcar + βtravel_time * CAR_TT / 100 + βcost * CAR_CO / 100

            αswissmetro = 0, fixed  # fix swissmetro ASC to zero
        end),
        :CHOICE,
        data,
        availability=[
            1 => :avtr,
            2 => :avsm,
            3 => :avcar,
        ]
    )
end

main()


text:
# Printing using `AbstractTrees.print_tree`
# I implemented this for Convex.jl and it seemed nice, so let's try it here too.

# A simple wrapper to stop `AbstractTrees.children` from recursing
struct NoChildren
    x::Any
end
Base.show(io::IO, n::NoChildren) = show(io, n.x)

AbstractTrees.children(q::NamedQuery) = (NoChildren(q.metadata), q.query)
AbstractTrees.children(q::NamedMatch) = (NoChildren(q.metadata), q.match)

AbstractTrees.children(q::QueryMatch) = (q.query, q.document)

AbstractTrees.children(q::Or) = q.subqueries
AbstractTrees.children(q::Query) = tuple()
AbstractTrees.children(q::FuzzyQuery) = tuple()

AbstractTrees.printnode(io::IO, q::Or) = print(io, "Or")
AbstractTrees.printnode(io::IO, q::NamedQuery) = print(io, "NamedQuery")
AbstractTrees.printnode(io::IO, q::NamedMatch) = print(io, "NamedMatch")

# Don't show all those types
function AbstractTrees.printnode(io::IO, q::FuzzyQuery)
    print(io, "FuzzyQuery(")
    show(io, q.text)
    print(io, ", ", q.dist, ", ", q.threshold, ")")
    return nothing
end

function print_tree_rstrip(io::IO, x)
    str = sprint(AbstractTrees.print_tree, x)
    print(io, rstrip(str))
    return nothing
end

"""
    Base.show(io::IO, q::Union{Or,NamedQuery,NamedMatch})

[`Or`]@(ref), [`NamedQuery`](@ref), and [`NamedMatch`](@ref) are
pretty-printed as trees.

## Example

```jldoctest
julia> q = Query("a") | Query("b")
Or
├─ Query("a")
└─ Query("b")

julia> named = NamedQuery(q, "is it a or is it b?")
NamedQuery
├─ (query_name = "is it a or is it b?",)
└─ Or
   ├─ Query("a")
   └─ Query("b")

julia> match(named, Document("It is a!"))
NamedMatch
├─ (query_name = "is it a or is it b?",)
└─ QueryMatch with distance 0 at indices 7:7.
   ├─ Query("a")
   └─ Document with text "It is a  ". Metadata: NamedTuple()

```
"""
Base.show(io::IO, q::Union{Or,NamedQuery,NamedMatch}) = print_tree_rstrip(io, q)

# `text_left_endpoint` and `text_right_endpoint` could be combined
# but it would probably be too clever... let's just duplicate it

function text_left_endpoint(text::AbstractString, start=1; approx_length=15)
    L = max(prevind(text, start, approx_length), firstindex(text))

    if L == firstindex(text)
        return L, ""
    end

    # Not at the start, so look for a word break
    space_inds = findprev(" ", text, L)

    if space_inds !== nothing
        L_space = nextind(text, last(space_inds))

        # make sure we don't have to go too far to get to the space
        if start - L_space <= 2 * approx_length
            L = L_space
        end
    end

    return L, "…"
end

function text_right_endpoint(text::AbstractString, start=1; approx_length=15)
    R = min(nextind(text, start, approx_length), lastindex(text))

    if R == lastindex(text)
        return R, ""
    end

    space_inds = findnext(" ", text, R)

    if space_inds !== nothing
        R_space = prevind(text, first(space_inds))
        if R_space - start <= 2 * approx_length
            R = R_space
        end
    end

    return R, "…"
end

"""
    Base.show(io::IO, D::Document)

Pretty-prints a [`Document`](@ref).

# Example

```jldoctest
julia> Document("Doc 1", (; doc_idx = 1)) # short documents print completely
Document with text "Doc 1 ". Metadata: (doc_idx = 1,)

julia> Document("This is a longer document! Lots of words here.", (; doc_idx = 1)) # longer documents print truncated
Document starting with "This is a longer…". Metadata: (doc_idx = 1,)

```
"""
function Base.show(io::IO, D::Document)
    n, rdots = text_right_endpoint(D.text)

    show_full_text = (length(D.text) < 15) || rdots == ""
    if show_full_text
        print(io, "Document with text ")
        show(io, D.text)
    else
        print(io, "Document starting with ")
        print(io, "\"", D.text[1:n], rdots, "\"")
    end
    print(io, ". Metadata: ")
    show(io, D.metadata)
    return nothing
end

"""
    Base.show(io::IO, C::Corpus)

Pretty-prints a [`Corpus`](@ref).

## Example

```jldoctest
julia> C1 = Corpus([Document("Doc 1", (; doc_idx = 1)), Document("Doc 2", (; doc_idx = 2))], (; name = "Lots of docs"))
Corpus with 2 documents, each with metadata keys: (:doc_idx,)
Corpus metadata: (name = "Lots of docs",)

julia> C2 = Corpus([Document("a")], (; a = 1));

julia> [C1, C2]
2-element Array{Corpus,1}:
 Corpus with 2 documents, each with metadata keys: (:doc_idx,)
 Corpus with 1 documents, each with metadata keys: ()

```
"""
function Base.show(io::IO, C::Corpus{T,TR}) where {T,TR}
    print(io, "Corpus with ", length(C.documents), " documents, each with metadata keys: ")
    print(io, _nt_names(TR))
    return nothing
end

function Base.show(io::IO, ::MIME"text/plain", C::Corpus{T,TR}) where {T,TR}
    show(io, C)
    print(io, "\nCorpus metadata: ", C.metadata)
    return nothing
end

function printcontext(io::IO, m::QueryMatch; context=40)
    text = m.document.text
    L, ldots = text_left_endpoint(text, first(m.indices); approx_length=context)
    R, rdots = text_right_endpoint(text, last(m.indices); approx_length=context)

    LL = prevind(text, first(m.indices))
    RR = nextind(text, last(m.indices))
    print(io, ldots)
    if LL > L
        print(io, replace(text[L:LL], "\n" => "\\n"))
    end
    printstyled(io, text[m.indices]; bold=true)
    if RR < R
        print(io, replace(text[RR:R], "\n" => "\\n"))
    end
    print(io, rdots)
    return nothing
end

# `raw` due to escaped quotes in `sprint(explain...)`
@doc raw"""
    explain([io=stdout], match; context=40)

Prints a human-readable explanation of the match
and its context in the document in which it was found.

## Example

```jldoctest
julia> document = Document("The crabeating macacue ate a crab.")
Document starting with "The crabeating macacue…". Metadata: NamedTuple()

julia> query = augment(FuzzyQuery("crab-eating macaque"))
Or
├─ FuzzyQuery("crab eating macaque", DamerauLevenshtein{Nothing}(nothing), 2)
├─ FuzzyQuery("crabeating macaque", DamerauLevenshtein{Nothing}(nothing), 2)
├─ FuzzyQuery("crab eatingmacaque", DamerauLevenshtein{Nothing}(nothing), 2)
└─ FuzzyQuery("crabeatingmacaque", DamerauLevenshtein{Nothing}(nothing), 2)

julia> m = match(query, document)
QueryMatch with distance 1 at indices 5:22.

julia> explain(m)
The query "crabeating macaque" matched the text "The crabeating macacue ate a crab  " with distance 1.

julia> explain(m; context=5) # tweak the amount of context printed
The query "crabeating macaque" matched the text "The crabeating macacue ate…" with distance 1.

julia> sprint(explain, m) # to get the explanation as a string
"The query \\"crabeating macaque\\" matched the text \\"The crabeating macacue ate a crab  \\" with distance 1.\n"

julia> explain(match(Query("crab"), document)) # exact queries print slightly differently
The query "crab" exactly matched the text "The crabeating macacue ate a crab  ".

julia> explain(match(NamedQuery(Query("crab"), "crab query"), document)) # `NamedQuery`s print the same as their underlying query
The query "crab" exactly matched the text "The crabeating macacue ate a crab  ".

```
"""
explain

function explain(io::IO, m::QueryMatch{<:FuzzyQuery}; context=40)
    print(io, "The query ")
    show(io, m.query.text)
    print(io, " matched the text \"")
    printcontext(io, m; context=context)
    print(io, "\" with distance ", m.distance, ".\n")
    return nothing
end

function explain(io::IO, m::QueryMatch{Query}; context=40)
    print(io, "The query ")
    show(io, m.query.text)
    print(io, " exactly matched the text \"")
    printcontext(io, m; context=context)
    print(io, "\".\n")
    return nothing
end

# in case we no longer have a `QueryMatch` but at least have a row-like
# object with the right fields.
function explain(io::IO, row; context=40)
    match = QueryMatch(row.query, row.document, row.distance, row.indices)
    return explain(io, match; context=context)
end

explain(m; context=40) = explain(stdout, m; context=context)

"""
    Base.show(io::IO, m::QueryMatch)

Pretty-prints a [`QueryMatch`](@ref).

## Example

```jldoctest
julia> match(Query("claws"), Document("Lemurs have nails instead of claws."))
QueryMatch with distance 0 at indices 30:34.

```
"""
function Base.show(io::IO, m::QueryMatch)
    return print(io, "QueryMatch with distance ", m.distance, " at indices ", m.indices,
                 ".")
end

max_stars_repo_name: beacon-biosignals/KeywordSearch

max_issues_repo_name: beacon-biosignals/KeywordSearch

max_forks_repo_name: beacon-biosignals/KeywordSearch


text:
using Printf
using ProgressMeter

push!(LOAD_PATH, @__DIR__)
using PCREO

function count_database_entries()
    result = 0
    dirnames = readdir(PCREO_DATABASE_DIRECTORY)
    prog = Progress(length(dirnames); desc="Scanning PCREO database... ")
    for dirname in dirnames
        result += length(readdir(joinpath(PCREO_DATABASE_DIRECTORY, dirname)))
        next!(prog)
    end
    println("Found $result database entries.")
    return result
end

function verify_unique_energies()
    count = count_database_entries()
    prog = Progress(count; desc="Verifying energy uniqueness... ")
    for dirname in readdir(PCREO_DATABASE_DIRECTORY)
        data = PCREORecord[]
        for fname in readdir(joinpath(PCREO_DATABASE_DIRECTORY, dirname))
            push!(data, PCREORecord(joinpath(
                PCREO_DATABASE_DIRECTORY, dirname, fname)))
            next!(prog)
        end
        energies = [record.energy for record in data]
        @assert length(unique!(energies)) == 1
    end
    println("Verified that all isometric database" *
            " entries have identical energies.")
end

function main()
    # verify_unique_energies()
    summary_data = Dict{Int,Vector{Tuple{Float64,String,Int}}}()
    @showprogress for dirname in readdir(PCREO_DATABASE_DIRECTORY)
        dirpath = joinpath(PCREO_DATABASE_DIRECTORY, dirname)
        filenames = readdir(dirpath)
        count = length(filenames)
        record = PCREORecord(joinpath(dirpath, first(filenames)))
        if record.num_points in keys(summary_data)
            push!(summary_data[record.num_points],
                  (record.energy, dirname, count))
        else
            summary_data[record.num_points] = [(record.energy, dirname, count)]
        end
    end
    for k in sort!(collect(keys(summary_data)))
        v = summary_data[k]
        sort!(v)
        total_count = sum(n for (_, _, n) in v)
        if length(v) > 1
            total_energy = sum(e * n for (e, _, n) in v)
            mean_energy = total_energy / total_count
            total_deviation = sum((e - mean_energy)^2 * n for (e, _, n) in v)
            stdev = sqrt(total_deviation / (total_count - 1))
        else
            mean_energy = v[1][1]
            stdev = 1.0
        end
        println("\n", k, " =======================================",
                "============================= ", total_count)
        for (e, s, n) in v
            z = (e - mean_energy) / stdev
            @printf("% 8.4f    %s    %4d\n", z, s, n)
        end
    end
end

main()


text:
using Base: tail
using ChainRules
using ChainRulesCore

"""
    ∂⃖{N}

∂⃖{N} is the reverse-mode AD optic functor of order `N`. A call
`(::∂⃖{N})(f, args...)` corresponds to  ∂⃖ⁿ f(args...) in the linear encoding
of an N-optic (see the terminology guide for definitions of these terms).

In general `(::∂⃖{N})(f, args...)` will return a tuple of the original primal
value `f(args...)` (in rare cases primitives may modify the primal value -
in general we will ignore this rare complication for the purposes of clear
documentation) and an optic continuation `λ`. The interpretation of this
continuation depends on the order of the functor:

For example, ∂⃖{1} computes first derivatives. In particular, for a function `f`,
`∂⃖{1}(f, args...)` will return the tuple `(f(args...), f⋆)` (read "f upper-star").

"""
struct ∂⃖{N}; end
const ∂⃖¹ = ∂⃖{1}()

(::Type{∂⃖})(args...) = ∂⃖¹(args...)

"""
    ∂☆{N}

∂☆{N} is the forward-mode AD functor of order `N`. A call
`(::∂☆{N})(f, args...)` evaluating a function `f: A -> B` is lifted to its
pushforward on the N-th order tangent bundle `f⋆: Tⁿ A -> Tⁿ B`.
"""
struct ∂☆{N}; end
const ∂☆¹ = ∂☆{1}()

"""
    dx(x)

dx represents the trival differential one-form of a one dimensional
Riemannian manifold `M`. In particular, it is a section of the cotangent bundle
of `M`, meaning it may be evaluted at a point `x` of `M` to obtain an element
of the cotangent space `T*ₓ M` to `M` at `x`. We impose no restrictions on
the representations of either the manifold itself or the cotangent space.

By default, the only implementation provided identifies T*ₓ ℝ ≃ ℝ, keeping
watever type is used to represent ℝ. i.e.
```julia
dx(x::Real) = one(x)
```

However, users may provide additional overloads for custom representations of
one dimensional Riemannian manifolds.
"""
dx(x::Real) = one(x)
dx(x::Complex) = error("Tried to take the gradient of a complex-valued function.")
dx(x) = error("Cotangent space not defined for `$(typeof(x))`. Try a real-valued function.")

"""
    ∂x(x)

For `x` in a one dimensional manifold, map x to the trivial, unital, 1st order
tangent bundle. It should hold that `∀x ⟨∂x(x), dx(x)⟩ = 1`
"""
∂x(x::Real) = TangentBundle{1}(x, (one(x),))
∂x(x) = error("Tangent space not defined for `$(typeof(x)).")

struct ∂xⁿ{N}; end

(::∂xⁿ{N})(x::Real) where {N} = TaylorBundle{N}(x, (one(x), (zero(x) for i = 1:(N-1))...,))
(::∂xⁿ)(x) = error("Tangent space not defined for `$(typeof(x)).")

function ChainRules.rrule(∂::∂xⁿ, x)
    ∂(x), Δ->(NoTangent(), Δ.primal)
end

"""
    ∇(f, args...)

Computes the gradient ∇f(x, y, z...) (at (x, y, z...)). In particular, the
return value will be a tuple of partial derivatives `(∂f/∂x, ∂f/∂y, ∂f/∂z...)`.

# Curried version

Alternatively, ∇ may be curried, essentially giving the gradient as a function:

## Examples

```jldoctest
julia> using Diffractor: ∇

julia> map(∇(*), (1,2,3), (4,5,6))
((4.0, 1.0), (5.0, 2.0), (6.0, 3.0))
```

# The derivative ∂f/∂f

Note that since in Julia, there is no distinction between functions and values,
there is in principle a partial derivative with respect to the function itself.
However, said partial derivative is dropped by this interface. It is however
available using the lower level ∂⃖ if desired. This interaction can also be used
to obtain gradients with respect to only some of the arguments by using a
closure:

∇((x,z)->f(x,y,z))(x, z) # returns (∂f/∂x, ∂f/∂z)

Though of course the same can be obtained by simply indexing the resulting
tuple (in well-inferred code there should not be a performance difference
between these two options).
"""
struct ∇{T}
    f::T
end

function Base.show(io::IO, f::∇)
    print(io, '∇')
    print(io, f.f)
end

function (f::∇)(args...)
    y, f☆ = ∂⃖(getfield(f, :f), args...)
    return tail(f☆(dx(y)))
end

# N.B: This means the gradient is not available for zero-arg function, but such
# a gradient would be guaranteed to be `()`, which is a bit of a useless thing
function (::Type{∇})(f, x1, args...)
    ∇(f)(x1, args...)
end

const gradient = ∇

# Star Trek has their prime directive. We have the...
abstract type AbstractPrimeDerivative{N, T}; end

# Backwards first order derivative
struct PrimeDerivativeBack{N, T} <: AbstractPrimeDerivative{N, T}
    f::T
end
Base.show(io::IO, f::PrimeDerivativeBack{N}) where {N} = print(io, f.f, "'"^N)

# This improves performance for nested derivatives by short cutting some
# recursion into the PrimeDerivative constructor
@Base.pure minus1(N) = N - 1
@Base.pure plus1(N) = N + 1
lower_pd(f::PrimeDerivativeBack{N,T}) where {N,T} = PrimeDerivativeBack{minus1(N),T}(getfield(f, :f))
lower_pd(f::PrimeDerivativeBack{1}) = getfield(f, :f)
raise_pd(f::PrimeDerivativeBack{N,T}) where {N,T} = PrimeDerivativeBack{plus1(N),T}(getfield(f, :f))

ChainRulesCore.rrule(::typeof(lower_pd), f) = lower_pd(f), Δ->(ZeroTangent(), Δ)
ChainRulesCore.rrule(::typeof(raise_pd), f) = raise_pd(f), Δ->(ZeroTangent(), Δ)

PrimeDerivativeBack(f) = PrimeDerivativeBack{1, typeof(f)}(f)
PrimeDerivativeBack(f::PrimeDerivativeBack{N, T}) where {N, T} = raise_pd(f)

function (f::PrimeDerivativeBack)(x)
    z = ∂⃖¹(lower_pd(f), x)
    y = getfield(z, 1)
    f☆ = getfield(z, 2)
    return getfield(f☆(dx(y)), 2)
end

# Forwards primal derivative
struct PrimeDerivativeFwd{N, T}
    f::T
end

PrimeDerivativeFwd(f) = PrimeDerivativeFwd{1, typeof(f)}(f)
PrimeDerivativeFwd(f::PrimeDerivativeFwd{N, T}) where {N, T} = raise_pd(f)

lower_pd(f::PrimeDerivativeFwd{N,T}) where {N,T} = (error(); PrimeDerivativeFwd{minus1(N),T}(getfield(f, :f)))
raise_pd(f::PrimeDerivativeFwd{N,T}) where {N,T} = PrimeDerivativeFwd{plus1(N),T}(getfield(f, :f))

(f::PrimeDerivativeFwd{0})(x) = getfield(f, :f)(x)

function (f::PrimeDerivativeFwd{1})(x)
    z = ∂☆¹(ZeroBundle{1}(getfield(f, :f)), ∂x(x))
    z.partials[1]
end

function (f::PrimeDerivativeFwd{N})(x) where N
    z = ∂☆{N}()(ZeroBundle{N}(getfield(f, :f)), ∂xⁿ{N}()(x))
    z[TaylorTangentIndex(N)]
end

# Polyalgorithm prime derivative
struct PrimeDerivative{N, T}
    f::T
end

function (f::PrimeDerivative{N, T})(x) where {N, T}
    # For now, this is backwards mode, since that's more fully implemented
    return PrimeDerivativeBack{N, T}(f.f)(x)
end

"""
    f'

This is a convenience syntax for taking the derivative of a function f: ℝ -> ℝ.
In particular, for such a function f'(x) will be the first derivative of `f`
at `x` (and similar for `f''(x)` and second derivatives and so on.)

Note that the syntax conflicts with the Base definition for the adjoint of a
matrix and thus is not enabled by default. To use it, add the following to the
top of your module:

```julia
using Diffractor: var"'"
```

It is also available using the @∂ macro:
```julia
@∂ f'(x)
```
"""
var"'"(f) = PrimeDerivativeBack(f)

"""
    @∂

Convenice macro for writing partial derivatives. E.g. The expression:

```julia
@∂ f(∂x, ∂y)
```

Will compute the partial derivative ∂^2 f/∂x∂y at `(x, y)``. And similarly

```julia
@∂ f(∂²x, ∂y)
```

will compute the derivative `∂^3 f/∂x^2 ∂y` at `(x,y)`.
"""
macro ∂(expr)
    error("Write me")
end


text:
  # Test routine to compare scaling for a number of SDP problems (partially badly scaled)
workspace()
include("../../../src/Solver.jl")

using OSSDP, Base.Test

rng = Random.MersenneTwister(12345)

nn = 20

resIter = zeros(20,2)
resCost = zeros(20,2)

for iii =1:1:nn
  # generate problem data
  m = 1
  r = 8
  n = r^2
  A = Symmetric(100*sprand(r,r,0.5))
  A = vec(A)'
  # create strictly feasible primal point
  X = Symmetric(10*sprand(r,r,0.5))
  X = X + (-minimum(eigs(X)[1])+1)*speye(r)
  b = A*vec(X)
  b = [b]

  # create strictly feasibile dual point
  y = rand(m,1)
  S = Symmetric(10*sprand(r,r,0.5))
  S = S + (-minimum(eigs(S)[1])+1)*speye(r)
  c = vec(S) + A'*y

  P = zeros(n,n)

  K = Cone(0,0,[],[n])

  # Solve with OSSDP
  settings = OSSDPSettings(rho=1.0,sigma=1.0,alpha=1.6,max_iter=2500,verbose=false,check_termination=1,scaling = 0,eps_abs = 1e-3,eps_rel=1e-4)
  resOSSDP_unscaled,nothing = OSSDP.solve(P,c,A,b,K,settings)
  settings = OSSDPSettings(rho=1.0,sigma=1.0,alpha=1.6,max_iter=2500,verbose=false,check_termination=1,scaling = 10,eps_abs = 1e-3,eps_rel=1e-4)
  resOSSDP_scaled,nothing = OSSDP.solve(P,c,A,b,K,settings)

  resCost[iii,:] = [resOSSDP_unscaled.cost resOSSDP_scaled.cost]
  resIter[iii,:] = [resOSSDP_unscaled.iter resOSSDP_scaled.iter]
  println("$(iii)/$(nn) completed!")
end

# @testset "Lasso Scaling" begin
#     @test abs(resOSSDP_unscaled.cost-resOSQP1_unscaled.info.obj_val) < 1e-2
#     @test norm(resOSSDP_unscaled.x[1:n]-resOSQP1_unscaled.x[1:n],Inf) < 1e-2
#     @test norm(resOSSDP_unscaled.x[1:n]-resOSQP2_unscaled.x[1:n],Inf) < 1e-2
#     @test norm(resOSQP1_unscaled.x[1:n]-resOSQP2_unscaled.x[1:n],Inf) < 1e-2
#   end


text:
"""
    single_phase_label(mixture, cond)

Li's method for single-phase labeling of a mixture. Estimate of pure vapor/liquid.

Returns a vapor fraction that is either 1.0 (=pure vapor) or 0.0 (=pure liquid).
"""
function single_phase_label(mixture, cond)
    # Li's method for phase labeling.
    z = cond.z
    T_c = 0.0
    V_c = 0.0
    for (i, m) in enumerate(mixture.properties)
        V = m.V_c*z[i]
        T_c += m.T_c*V
        V_c += V
    end
    return Float64(cond.T < T_c)
end

"""
    lbc_viscosity(eos, p, T, ph; <keyword arguments>)

Compute the viscosity of a mixture using the [Lohrenz-Bray-Clark](https://doi.org/10.2118/915-PA) correlation.
"""
function lbc_viscosity(eos, p, T, ph; coeff = (0.1023, 0.023364, 0.058533, -0.040758, 0.0093324), shift = -1e-4)
    molfactor = 1000;
    Rankine = 5/9
    psia = 6.894757293168360e+03

    rho = 1/molar_volume(eos, p, T, ph)
    z = ph.mole_fractions

    properties = eos.mixture.properties
    # for i = 1:ncomp
    P_pc = 0.0
    T_pc = 0.0
    Vc = 0.0
    mwc = 0.0

    a = 0.0
    b = 0.0
    for i in eachindex(properties)
        prop = properties[i]
        zi = z[i]
        mw = prop.mw
        T_c = prop.T_c
        p_c = prop.p_c
        V_c = prop.V_c

        mwi = sqrt(molfactor*mw)

        mwc += zi*mw
        P_pc += zi*p_c
        T_pc += zi*T_c
        Vc += zi*V_c

        tr = T/T_c
        Tc = T_c/Rankine;
        Pc = p_c/psia;

        e_i = (5.4402*Tc.^(1/6))./(mwi.*Pc.^(2/3).*(1e-3));

        large = tr > 1.5;
        if large
            mu_i = (17.78e-5.*(4.58*tr - 1.67).^0.625)./e_i;
        else
            mu_i = 34e-5.*tr.^(0.94)./e_i;
        end

        a = a + zi.*mu_i.*mwi;
        b = b + zi.*mwi;
    end
    mu_atm = a./b;
    e_mix = 5.4402*(T_pc./Rankine).^(1/6)./((molfactor*mwc).^(1/2).*(P_pc./psia).^(2/3).*(1e-3));
    rhor = Vc.*rho;
    # From Jossi et al
    # coeffs = [0.1023, 0.023364, 0.058533, -0.040758, 0.0093724, 1e-4];
    # From LBC paper
    # coeffs = [0.1023, 0.023364, 0.058533, -0.040758, 0.0093324, -1e-4];
    corr = 0.0
    for (i, c) in enumerate(coeff)
        corr += c*rhor^(i-1)
    end
    mu = mu_atm + (corr^4 + shift)./e_mix;
    return mu
end

function molar_volume(R, p, T, Z)
    return R*T*Z/p
end


text:
#!/usr/bin/env julia

#load path to qjulia home directory
push!(LOAD_PATH, joinpath(@__DIR__, "..", "matrix"))
push!(LOAD_PATH, joinpath(@__DIR__, "..", "core"))

using MatrixBase
using MatrixMarket
using LinearAlgebra
using Random

using QJuliaSolvers
using QJuliaEnums
using MPI
using Printf

#matrix_path = "/home/astrel/data/bcsstk15/bcsstk15.mtx"
matrix_path = "/home/alex/data/bcsstk03/bcsstk03.mtx"
#matrix_path = "/home/astrel/data/nasa2146/nasa2146.mtx"
#matrix_path = "/home/astrel/data/nasa2910/nasa2910.mtx"
#matrix_path = "/home/astrel/data/nasa4704/nasa4704.mtx"
#matrix_path = "/home/astrel/data/smt/smt.mtx"
#matrix_path = "/home/alex/data/nos4/nos4.mtx"

#initialize MPI
MPI.Init()

Random.seed!(2018)

prec        = Float64
prec_sloppy = Float64

setprecision(BigFloat, 128) #128 bit

csrM       = MatrixBase.CSRMat{prec}(matrix_path)
#MatrixBase.rescale(csrM, 0.5)
csrMsloppy = prec_sloppy != prec ? MatrixBase.CSRMat{prec_sloppy}(matrix_path) : csrM

MatrixBase.print_CSRMat_info(csrM)

data_type = typeof(csrM.csrVals[1])

b = Vector{data_type}(undef, csrM.Dims[2])
x = Vector{data_type}(undef, csrM.Dims[2])
z = Vector{data_type}(undef, csrM.Dims[2])

x .= 0.0
z .= 0.0
#b .= rand(data_type, csrM.Dims[2])

M(b,x)        = MatrixBase.csrmv(b, csrM, x)

if prec_sloppy != prec
  Msloppy(b,x)  = MatrixBase.csrmv(b, csrMsloppy, x)
else
  Msloppy = M
end
#or just:
#M=MatrixMarket.mmread(matrix_path)
# check matrix vector operation:
println("Begin check")
w64 = ones(data_type, csrM.Dims[2])
v64 = ones(data_type, csrM.Dims[2])
data_type_sloppy = typeof(csrMsloppy.csrVals[1])
w32 = ones(data_type_sloppy, csrM.Dims[2])
v32 = ones(data_type_sloppy, csrM.Dims[2])

M(v64, w64)
w64 .=@. abs.(v64)
Ainfnorm = findmax(w64)[1]
println("\nExtimated matrix inf norm: \n", Ainfnorm)

ϵ64 = eps(data_type)
ϵ32 = eps(data_type_sloppy)

println("Maximum error bounds: ", Ainfnorm*ϵ64, " for FP and ", Ainfnorm*ϵ32, " for SP")
println("Check unit sources.")
w64 .=@. 1.0
v64 .=@. 0.0
M(v64, w64)
w32 .=@. 1.0
v32 .=@. 0.0
Msloppy(v32, w32)

w64 .=@. v64 - v32
@printf("Difference norm %1.10le,\t(%1.10le vs %1.10le)\n", norm(w64), norm(v32), norm(v64))

println("Check random sources.")

w64 .= rand(data_type, csrM.Dims[2])
v64 .=@. 0.0
M(v64, w64)
w32 .=@. w64
v32 .=@. 0.0
Msloppy(v32, w32)

w64 .=@. v64 - v32
@printf("Difference norm %1.10le,\t(%1.10le vs %1.10le)\n", norm(w64), norm(v32), norm(v64))

println("Check random source cast.")

w64 .= rand(data_type, csrM.Dims[2])
v64 .=@. 0.0
w32 .=@. w64
v32 .=@. 0.0
Msloppy(v32, w32)
w64 .=@. w32
M(v64, w64)

w64 .=@. v64 - v32
@printf("Difference norm %1.10le,\t(%1.10le vs %1.10le)\n", norm(w64), norm(v32), norm(v64))

println("Finish check \n")

solv_param = QJuliaSolvers.QJuliaSolverParam_qj()
# Set up parameters
solv_param.dtype                  = prec
solv_param.dtype_sloppy           = prec_sloppy
solv_param.inv_type               = QJuliaEnums.QJULIA_PIPEPCG_INVERTER
#solv_param.inv_type               = QJuliaEnums.QJULIA_PCG_INVERTER
solv_param.inv_type_precondition  = QJuliaEnums.QJULIA_INVALID_INVERTER
solv_param.tol                    = 1e-10
solv_param.delta                  = 1e+1
#
solv_param.maxiter                = 1000
solv_param.Nsteps                 = 2

do_unit_source = true
if do_unit_source == true
  z .=@. 1.0 / sqrt(csrM.Dims[2])
else
  z .= rand(data_type, csrM.Dims[2])
end

M(b,z)

println("Initial residual:  ", norm(b))

@time QJuliaSolvers.solve(x, b, M, Msloppy, solv_param)

#compute error:
z .=@. z - x
println("Solution error:  ", norm(z))

#compute residual:
M(z, x)
z .= b - z

println("Solution norm: ", norm(x), ", absolute residual: ", norm(z), " relative ", norm(z) / norm(b))

MPI.Finalize()

max_stars_repo_name: alexstrel/QJulia

max_issues_repo_name: alexstrel/QJulia

max_forks_repo_name: alexstrel/QJulia


text:
#!/usr/bin/env julia
# This file is a part of Julia. License is MIT: http://julialang.org/license

const options =
[
    "--cflags",
    "--ldflags",
    "--ldlibs"
];

threadingOn() = ccall(:jl_threading_enabled, Cint, ()) != 0

function imagePath()
    opts = Base.JLOptions()
    unsafe_string(opts.image_file)
end

function libDir()
    return if ccall(:jl_is_debugbuild, Cint, ()) != 0
        dirname(abspath(Libdl.dlpath("libjulia-debug")))
    else
        dirname(abspath(Libdl.dlpath("libjulia")))
    end
end

function includeDir()
    joinpath(match(r"(.*)(bin)",JULIA_HOME).captures[1],"include","julia")
end

function unixInitDir()
    filePart = split(imagePath(),"/")[end]
    return match(Regex("(.*)(/julia/$filePart)"),imagePath()).captures[1]
end

function windowsInitDir()
    if imagePath()[end-1:end] == "ji"
        return match(r"(.*)(\\julia\\sys.ji)",imagePath()).captures[1]
    else
        return match(r"(.*)(\\julia\\sys.dll)",imagePath()).captures[1]
    end
end

function initDir()
    if is_unix()
        return unixInitDir()
    else
        return windowsInitDir()
    end
end

function ldflags()
    replace("""-L$(libDir())""","\\","\\\\")
end

function ldlibs()
    libname = if ccall(:jl_is_debugbuild, Cint, ()) != 0
        "julia-debug"
    else
        "julia"
    end
    if is_unix()
        return replace("""-Wl,-rpath,$(libDir()) -l$libname""","\\","\\\\")
    else
        return "-l$libname -lopenlibm"
    end
end

function cflags()
    arg1 = replace(initDir(),"\\","\\\\\\\\")
    arg2 = replace(includeDir(),"\\","\\\\")
    threading_def = threadingOn() ? "-DJULIA_ENABLE_THREADING=1 " : ""
    if is_unix()
        return """$(threading_def)-fPIC -DJULIA_INIT_DIR=\\"$arg1\\" -I$arg2"""
    else
        return """$(threading_def)-DJULIA_INIT_DIR=\\"$arg1\\" -I$arg2"""
    end
end

function check_args(args)
    checked = intersect(args,options)
    if length(checked) == 0 || length(checked) != length(args)
        println(STDERR,"Usage: julia-config [",reduce((x,y)->"$x|$y",options),"]")
        exit(1)
    end
end

function main()
    check_args(ARGS)
    for args in ARGS
        if args == "--ldflags"
            println(ldflags())
        elseif args == "--cflags"
            println(cflags())
        elseif args == "--ldlibs"
            println(ldlibs())
        end
    end
end

main()


text:
using FeatherLib, Missings, Dates, CategoricalArrays, Random, Arrow
using Test

temps = []

@testset "FeatherLib" begin

include("test_readwrite.jl")
include("test_arrow.jl")

GC.gc(); GC.gc()
for t in temps
    try
        rm(t)
    catch
        GC.gc()
        try
            rm(t)
        catch
        end
    end
end

# issue #34
# data = DataFrame(A=Union{Missing, String}[randstring(10) for i ∈ 1:100], B=rand(100))
# data[2, :A] = missing
# Feather.write("testfile.feather", data)
# dfo = Feather.read("testfile.feather")
# @test size(Data.schema(dfo)) == (100, 2)
# GC.gc();
# rm("testfile.feather")

# @testset "PythonRoundtrip" begin
# try
#     println("Generate a test.feather file from python...")
#     run(`docker cp runtests.py feathertest:/home/runtests.py`)
#     run(`docker exec feathertest python /home/runtests.py`)

#     println("Read test.feather into julia...")
#     run(`docker cp feathertest:/home/test.feather test.feather`)
#     df = Feather.read("test.feather")

#     dts = [Dates.DateTime(2016,1,1), Dates.DateTime(2016,1,2), Dates.DateTime(2016,1,3)]

#     @test df[:Autf8][:] == ["hey","there","sailor"]
#     @test df[:Abool][:] == [true, true, false]
#     @test df[:Acat][:] == categorical(["a","b","c"])  # these violate Arrow standard by using Int8!!
#     @test df[:Acatordered][:] == categorical(["d","e","f"])  # these violate Arrow standard by using Int8!!
#     @test convert(Vector{Dates.DateTime}, df[:Adatetime][:]) == dts
#     @test isequal(df[:Afloat32][:], [1.0, missing, 0.0])
#     @test df[:Afloat64][:] == [Inf,1.0,0.0]

#     df_ = Feather.read("test.feather"; use_mmap=false)

#     println("Writing test2.feather from julia...")
#     Feather.write("test2.feather", df)
#     df2 = Feather.read("test2.feather")

#     @test df2[:Autf8][:] == ["hey","there","sailor"]
#     @test df2[:Abool][:] == [true, true, false]
#     @test df2[:Acat][:] == categorical(["a","b","c"])  # these violate Arrow standard by using Int8!!
#     @test df2[:Acatordered][:] == categorical(["d","e","f"])  # these violate Arrow standard by using Int8!!
#     @test convert(Vector{Dates.DateTime}, df2[:Adatetime][:]) == dts
#     @test isequal(df2[:Afloat32][:], [1.0, missing, 0.0])
#     @test df2[:Afloat64][:] == [Inf,1.0,0.0]

#     println("Read test2.feather into python...")
#     @test (run(`docker cp test2.feather feathertest:/home/test2.feather`); true)
#     @test (run(`docker cp runtests2.py feathertest:/home/runtests2.py`); true)
#     @test (run(`docker exec feathertest python /home/runtests2.py`); true)
# finally
#     run(`docker stop feathertest`)
#     run(`docker rm feathertest`)
#     rm("test.feather")
#     rm("test2.feather")
# end

# end

end


text:
function step_forward!(m::ClimateModel, Δt::Float64)
    m.domain.present_year = m.domain.present_year + Δt
    m.name = string(Int64(round(m.domain.present_year)))
end

function add_emissions_bump!(m::ClimateModel, Δt::Float64, Δq::Float64; present_year = m.domain.present_year)

    present_idx = argmin(abs.(t(m) .- (present_year .+ Δt)))

    future = (t(m) .>= present_year)
    near_future = future .& (t(m) .<= present_year + Δt)
    near_future1 = near_future .& (t(m) .< present_year + Δt/2)
    near_future2 = near_future .& (t(m) .>= present_year + Δt/2)

    new_emissions = m.economics.baseline_emissions
    new_emissions[near_future1] .+= (
        Δq * (t(m) .- present_year) / (Δt/2.)
    )[near_future1]
    new_emissions[near_future2] .+= (
        Δq * (1. .- (t(m) .- (present_year .+ Δt/2.)) / (Δt/2.))
    )[near_future2]

    m.economics.baseline_emissions = new_emissions
end


text:
"""
Complex numbers as 2-vectors, enough for polynomials.
"""
struct Compl{T} <: Number
    re::T
    im::T
end

Base.promote_rule{T, S<:Real}(::Type{Compl{T}}, ::Type{S}) = Compl{T}
Base.convert{T}(::Type{Compl{T}}, x::Real) = Compl{T}(x, 0)

Base.show(io::IO, c::Compl) = println(io, c.re, " + ", c.im, "im")

Base.:+{T}(a::Compl{T}, b::Compl{T}) = Compl{T}(a.re+b.re, a.im+b.im)
Base.:-{T}(a::Compl{T}, b::Compl{T}) = Compl{T}(a.re-b.re, a.im-b.im)

Base.:*{T}(a::Compl{T}, b::Compl{T}) = Compl{T}(a.re*b.re - a.im*b.im, a.re*b.im + a.im*b.re)

Base.:*{T}(α::Number, z::Compl{T}) = Compl{T}(α*z.re, α*z.im)

Base.one{T}(z::Compl{T}) = Compl{T}(one(T), zero(T))

Base.copy{T}(z::Compl{T}) = Compl{T}(z.re, z.im)

"""
Takes a complex (polynomial) function f and returns a function g:R^2 -> R^2
that implements it.
"""
function realify(f)

    function g(x)
        z = Compl(x[1], x[2])
        z2 = f(z)
        SVector(z2.re, z2.im)
    end

    return g

end