text:
############
# Starring #
############

#= Seems like there is no equivalent of this method to get the watchers of repos
function stargazers(repo; options...)
    ## results, page_data = gh_get_paged_json("/repos/$(name(repo))/stargazers"; options...)
    results, page_data = gh_get_paged_json("/api/v3/projects/$(get(repo.id))/star"; options...)
    return map(Owner, results), page_data
end
=#

function starred(user; options...)
    ## results, page_data = gh_get_paged_json("/api/v3/users/$(name(user))/starred"; options...)
    results, page_data = gh_get_paged_json("/api/v3/projects/starred"; options...)
    return map(Repo, results), page_data
end

## star(repo; options...) = gh_put("/user/starred/$(name(repo))"; options...)
star(repo; options...) = gh_post("/api/v3/projects/$(get(repo.id))/star"; options...)

## unstar(repo; options...) = gh_delete("/user/starred/$(name(repo))"; options...)
unstar(repo; options...) = gh_delete("/api/v3/projects/$(get(repo.id))/star"; options...)

############
# Watching #
############

#= Seems like there is no equivalent of these methods
function watchers(repo; options...)
    ## results, page_data = gh_get_paged_json("/repos/$(name(repo))/subscribers"; options...)
    results, page_data = gh_get_paged_json("/api/v3/projects/$(get(repo.id))/subscribers"; options...)
    return map(Owner, results), page_data
end

function watched(owner; options...)
    results, page_data = gh_get_paged_json("/api/v3/users/$(name(owner))/subscriptions"; options...)
    return map(Repo, results), page_data
end

## watch(repo; options...) = gh_put("/repos/$(name(repo))/subscription"; options...)
watch(repo; options...) = gh_put("/api/v3/projects/$(get(repo.id))/subscription"; options...)

## unwatch(repo; options...) = gh_delete("/repos/$(name(repo))/subscription"; options...)
unwatch(repo; options...) = gh_delete("/api/v3/projects/$(get(repo.id))/subscription"; options...)
=#


text:
# Intersect
# =========

immutable BAMIntersectionIterator{T}
    # data reader
    reader::BAMReader{T}

    # query interval
    refname::AbstractString
    interval::UnitRange{Int}
end

function Base.iteratorsize(::BAMIntersectionIterator)
    return Base.SizeUnknown()
end

function Base.eltype{T}(::Type{BAMIntersectionIterator{T}})
    return BAMRecord
end

function Base.show(io::IO, iter::BAMIntersectionIterator)
    print(
        io,
        summary(iter),
        "(<$(iter.refname):$(first(iter.interval))-$(last(iter.interval))>)")
end

type BAMIntersectionIteratorState
    # reference index
    refindex::Int

    # possibly overlapping chunks
    chunks::Vector{Chunk}

    # current chunk index
    chunkid::Int

    # pre-allocated record
    record::BAMRecord
end

function Base.start(iter::BAMIntersectionIterator)
    refindex = findfirst(iter.reader.refseqnames, iter.refname)
    if refindex == 0
        throw(ArgumentError("sequence name $(iter.refname) is not found in the header"))
    end
    @assert !isnull(iter.reader.index)
    chunks = overlapchunks(get(iter.reader.index).index, refindex, iter.interval)
    if !isempty(chunks)
        seek(iter.reader, first(chunks).start)
    end
    return BAMIntersectionIteratorState(refindex, chunks, 1, BAMRecord())
end

function Base.done(iter::BAMIntersectionIterator, state)
    while state.chunkid ≤ endof(state.chunks)
        chunk = state.chunks[state.chunkid]
        while virtualoffset(iter.reader.stream) < chunk.stop
            read!(iter.reader, state.record)
            if isoverlapping(state.record, state.refindex, iter.interval)
                return false
            end
        end
        state.chunkid += 1
        if state.chunkid ≤ endof(state.chunks)
            seek(iter.reader, state.chunks[state.chunkid].start)
        end
    end
    return true
end

function Base.next(iter::BAMIntersectionIterator, state)
    return copy(state.record), state
end

function Bio.Intervals.isoverlapping(
        rec::BAMRecord,
        refindex_::Integer,
        interval::UnitRange)
    return refindex(rec) == refindex_ &&
        leftposition(rec) ≤ last(interval) &&
        rightposition(rec) ≥ first(interval)
end

function Base.intersect(reader::BAMReader, interval::Interval)
    return intersect(reader, interval.seqname, interval.first:interval.last)
end

function Base.intersect(
        reader::BAMReader,
        refname::AbstractString,
        interval::UnitRange)
    return BAMIntersectionIterator(reader, refname, interval)
end

ext: jl


text:
"""
This module implements computation of kinship and other identity coefficients.
The kinship matrix is computed recursively by the standard formula.
The delta7 matrix is derived from the kinship matrix for non-inbred pedigrees.
The remaining identity coefficients of Jacquard are computed by simulation.
"""
module MendelKinship
#
# Required OpenMendel packages and modules.
#
using MendelBase
using SnpArrays
#
# Required external modules.
#
using PlotlyJS
using StatsBase
using DataFrames                        # From package DataFrames.
##
using CSV
##

export Kinship

"""
This is the wrapper function for the Kinship analysis option.
"""
function Kinship(control_file = ""; args...)
#
  const KINSHIP_VERSION :: VersionNumber = v"0.2.0"
#
# Print the logo. Store the initial directory.
#
  print(" \n \n")
  println("     Welcome to OpenMendel's")
  println("     Kinship analysis option")
  println("        version ", KINSHIP_VERSION)
  print(" \n \n")
  println("Reading the data.\n")
  initial_directory = pwd()
#
# The user specifies the analysis option via a set of keywords.
# Start the keywords at their default values.
#
  keyword = set_keyword_defaults!(Dict{AbstractString, Any}())
#
# Keywords unique to this analysis should be first defined here.
# Set their default values using the format:
# keyword["some_keyword_name"] = default_value
#
  keyword["repetitions"] = 1
  keyword["xlinked_analysis"] = false
  keyword["kinship_output_file"] = "kinship_file_output.txt"
  keyword["compare_kinships"] = false
  keyword["maf_threshold"] = 0.01
  keyword["grm_method"] = "MoM" # MoM is less sensitive to rare snps
  keyword["deviant_pairs"] = 0
  keyword["kinship_plot"] = ""
  keyword["z_score_plot"] = ""
#
# Process the run-time user-specified keywords that will control the
# analysis. This will also initialize the random number generator.
#
  process_keywords!(keyword, control_file, args)
#
# Check that the correct analysis option was specified.
#
  lc_analysis_option = lowercase(keyword["analysis_option"])
  if (lc_analysis_option != "" &&
      lc_analysis_option != "kinship")
     throw(ArgumentError(
       "An incorrect analysis option was specified.\n \n"))
  end
  keyword["analysis_option"] = "Kinship"
#
# Read the genetic data from the external files named in the keywords.
#
  (pedigree, person, nuclear_family, locus, snpdata, locus_frame,
    phenotype_frame, pedigree_frame, snp_definition_frame) =
    read_external_data_files(keyword)
#
# Execute the specified analysis.
#
  println(" \nAnalyzing the data.\n")
  if keyword["compare_kinships"]
    @assert snpdata.snps != 0 "Need Snp data to compare kinships"
    @assert snpdata.snps == length(unique(snpdata.snpid)) "non-unique snp ids"
    @assert person.people == length(unique(person.name)) "non-unique person names"
#
# Find each person's order on data entry before MendelBase permutes them.
#
    entryorder = convert(Vector{Int64}, pedigree_frame[:EntryOrder])
#
# Calculate empiric and theoretical kinships and append fisher's z-scores
# of these to the kinship frame.
#
    kinship_frame = compare_kinships(pedigree, person, snpdata, keyword, entryorder)
    kinship_frame = compute_fishers_z(kinship_frame)
#
# Create and save plots based on user's requests.
#
    name = make_plot_name(kinship_frame)
    if keyword["kinship_plot"] != ""
      my_compare_plot = make_compare_plot(kinship_frame, name)
      PlotlyJS.savefig(my_compare_plot, keyword["kinship_plot"] * ".html")
      println("Kinship plot saved.")
    end
    if keyword["z_score_plot"] != ""
      my_fisher_plot = plot_fisher_z(kinship_frame, name)
      PlotlyJS.savefig(my_fisher_plot, keyword["z_score_plot"] * ".html")
      println("Fisher's plot saved.")
    end
  else
#
# Otherwise, simply compute theoretical coefficients.
#
    kinship_frame = kinship_option(pedigree, person, keyword)
  end
#
# Save the kinship frame and write it to the designated file and REPL.
#
  name = convert(String, keyword["kinship_output_file"])
  CSV.write(name, kinship_frame)
  display(kinship_frame)
#
# Finish up by closing, and thus flushing, the output file.
# Return to the initial directory.
#
  println(" \n \nMendel's analysis is finished.\n")
  close(keyword["output_unit"])
  cd(initial_directory)
  return nothing
end # function Kinship

"""Compares theoretical and empirical kinship coefficients."""

function compare_kinships(pedigree::Pedigree, person::Person,
  snpdata::SnpData, keyword::Dict{AbstractString, Any}, entryorder::Vector{Int})
#
  people = person.people
  pedigrees = pedigree.pedigrees
  xlinked = keyword["xlinked_analysis"]
  maf_threshold = keyword["maf_threshold"]
  deviant_pairs = keyword["deviant_pairs"]
  if deviant_pairs == 0
    deviant_pairs = people^2
  end
#
# Create a vector of kinship matrices, one for each pedigree.
#
  kinship = Vector{Matrix{Float64}}(pedigrees)
##
## Read the full pedigree, including people with no snp information
##
  # ped_full = CSV.read("Ped10cextra.csv")
#
# Create a copy of snpdata.personid that orders people similarly to person.name.
# People are permuted so the parents come before children. This facilitates
# computation of theoretical kinships. In SnpArrays people are not permuted.
#
  personid_in_snpdata = copy(snpdata.personid)
  ipermute!(personid_in_snpdata, entryorder)
#
# Since people are reordered, compute identification maps to match
# positions of ids in two vectors of ids.
#
  name_to_id = indexin(person.name, personid_in_snpdata)
  id_to_name = indexin(personid_in_snpdata, person.name)
#
# Compute and transform the genetic relationship matrix.
#
  if keyword["grm_method"] == "GRM"
    GRM = grm(snpdata.snpmatrix, method=:GRM, maf_threshold=maf_threshold)
  else
    GRM = grm(snpdata.snpmatrix, method=:MoM)
  end
  clamp!(GRM, -0.99999, 0.99999)
  map!(x -> atanh(x), GRM, GRM) #fisher's transformation
#
# Loop over all pedigrees.
#
  for ped = 1:pedigrees
#
# Compute theoretical coefficients and adjust the GRM matrix.
#
    kinship[ped] = kinship_matrix(pedigree, person, ped, xlinked)
    q = pedigree.start[ped] - 1
    for j = 1:pedigree.individuals[ped]
      jj = name_to_id[j + q]
      if jj == 0; continue; end
      for i = j:pedigree.individuals[ped]
        ii = name_to_id[i + q]
        if ii == 0; continue; end
        GRM[ii, jj] = GRM[ii, jj] - atanh(kinship[ped][i, j])
        GRM[jj, ii] = GRM[ii, jj]
      end
    end
  end
#
# Find the indices of the most deviant pairs.
#
  p = selectperm(vec(GRM), 1:deviant_pairs, by = abs, rev = true)
  (rowindex, columnindex) = ind2sub((people, people), p)
#
# Enter the most deviant pairs in a data frame.
#
  kinship_frame = DataFrame(Pedigree1 = String[], Pedigree2 = String[],
    Person1 = String[], Person2 = String[], theoretical_kinship= Float64[],
    empiric_kinship=Float64[])
  r = 0.0
  for k = 1:deviant_pairs
    (ii, jj) = (rowindex[k], columnindex[k]) # index of kth largest deviation
    (i, j) = (id_to_name[ii], id_to_name[jj]) # maps snpid back to person name
    if j > i; continue; end
    (pedi, pedj) = (person.pedigree[i], person.pedigree[j]) # pedigrees of i and j
    if pedi == pedj
      q = pedigree.start[pedi] - 1
      r = GRM[ii, jj] + atanh(kinship[pedi][i - q, j - q])
      push!(kinship_frame, [pedigree.name[pedi], pedigree.name[pedj], person.name[i],
        person.name[j], kinship[pedi][i - q, j - q], tanh(r)])
    else
      push!(kinship_frame, [pedigree.name[pedi], pedigree.name[pedj], person.name[i],
        person.name[j], 0.0, tanh(GRM[ii, jj])])
    end
  end
  return kinship_frame
end

"""
This function orchestrates the computation via gene dropping
of the Kinship and Delta7 coefficients deterministically
and Jacquard's 9 identity coefficients stochastically.
The results are placed in a dataframe.
"""
function kinship_option(pedigree::Pedigree, person::Person,
  keyword::Dict{AbstractString, Any})

  pedigrees = pedigree.pedigrees
  repetitions = keyword["repetitions"]
  xlinked = keyword["xlinked_analysis"]
  #
  # Define dataframes for theoretical Kinship, Delta7, and Jacquard
  # coefficients.
  #
    # kinship_frame = DataFrame(Pedigree1 = String[], Pedigree2 = String[],
    # Person1 = String[], Person2 = String[], theoretical_kinship= Float64[],
    # empiric_kinship=Float64[])
  kinship_dataframe = DataFrame(
    ped1 = Int[], Pedigree = String[], Person1 = String[], Person2 = String[], Kinship = Float64[])
  if !xlinked
    delta_dataframe = DataFrame(ped2 = Int[], Pedigree = String[], Person1 = String[],
      Person2 = String[], Delta7 = Float64[])
  end
  coefficient_dataframe = DataFrame(
    ped3 = Int[], Pedigree = String[], Person1 = String[], Person2 = String[],
    delta1 = Float64[], delta2 = Float64[], delta3 = Float64[],
    delta4 = Float64[], delta5 = Float64[], delta6 = Float64[],
    delta7 = Float64[], delta8 = Float64[], delta9 = Float64[])
  #
  # Loop over all pedigrees.
  #
  for ped = 1:pedigrees
    #
    # Retrieve the various coefficients.
    #
    kinship = kinship_matrix(pedigree, person, ped, xlinked)
    if !xlinked
      delta7 = delta7_matrix(pedigree, person, kinship, ped)
    end
    coefficient = jacquard_coefficients(pedigree, person, ped,
      repetitions, xlinked)
    #
    # Insert the computed coefficients into the appropriate dataframes.
    #
    q = pedigree.start[ped] - 1
    for i = 1:pedigree.individuals[ped]
      for j = i:pedigree.individuals[ped]
        push!(kinship_dataframe, [ped, pedigree.name[ped],
          person.name[i + q], person.name[j + q], kinship[i, j]])
        if !xlinked
          push!(delta_dataframe, [ped, pedigree.name[ped],
            person.name[i + q], person.name[j + q], delta7[i, j]])
        end
        push!(coefficient_dataframe, [ped, pedigree.name[ped],
          person.name[i + q], person.name[j + q],
          coefficient[1, i, j], coefficient[2, i, j], coefficient[3, i, j],
          coefficient[4, i, j], coefficient[5, i, j], coefficient[6, i, j],
          coefficient[7, i, j], coefficient[8, i, j], coefficient[9, i, j]])
      end
    end
  end
  #
  # Join the separate dataframes, and sort the combined dataframe so
  # that relative pairs are lumped by pedigree.
  #
  if xlinked
    combined_dataframe = join(kinship_dataframe, coefficient_dataframe,
      on = [:Pedigree, :Person1, :Person2])
    sort!(combined_dataframe, cols = [:ped1, :Person1, :Person2])
    delete!(combined_dataframe, [:ped1, :ped3])
  else
    temp_dataframe = join(kinship_dataframe, delta_dataframe,
      on = [:Pedigree, :Person1, :Person2])
    combined_dataframe = join(temp_dataframe, coefficient_dataframe,
      on = [:Pedigree, :Person1, :Person2])
    sort!(combined_dataframe, cols = [:ped1, :Person1, :Person2])
    delete!(combined_dataframe, [:ped1, :ped2, :ped3])
  end
  #
  # Combined coefficient frame to a file and return.
  #
  return combined_dataframe
end # function kinship_option

"""
This function recursively computes a kinship matrix for each pedigree.
Children must appear after parents. Each pedigree occurs in a block
demarked by its start and finish and is defined implicitly by
the provided fathers and mothers. Person i is a pedigree founder
if and only if mother[i] = 0. X-linked inheritance is allowed.
"""
function kinship_matrix(pedigree::Pedigree, person::Person,
  ped::Int, xlinked::Bool)
  #
  # Allocate a kinship matrix and an offset.
  #
  start = pedigree.start[ped]
  finish = pedigree.finish[ped]
  kinship = zeros(pedigree.individuals[ped], pedigree.individuals[ped])
  q = start - 1
  for i = start:finish
    #
    # Initialize the kinship coefficients for founders.
    #
    j = person.mother[i]
    if j == 0
      if xlinked && person.male[i]
        kinship[i - q, i - q] = 1.0
      else
        kinship[i - q, i - q] = 0.5
      end
    else
      #
      # Compute the kinship coefficients of descendants by averaging.
      #
      if xlinked && person.male[i]
        kinship[i - q, i - q] = 1.0
        for m = start:i - 1
          kinship[i - q, m - q] = kinship[j - q, m - q]
          kinship[m - q, i - q] = kinship[i - q, m - q]
        end
      else
        k = person.father[i]
        kinship[i - q, i - q] = 0.5 + 0.5 * kinship[j - q, k - q]
        for m = start:i - 1
          kinship[i - q, m - q] = 0.5 * (kinship[j - q, m - q]
             + kinship[k - q, m - q])
          kinship[m - q, i - q] = kinship[i - q, m - q]
        end
      end
    end
  end
  #
  # Extend the kinship matrix to co-twins.
  #
  cotwin_extension!(kinship, pedigree, person, ped)
  return kinship
end # function kinship_matrix

"""
This function computes the Delta7 matrix for non-inbred pedigrees.
Children must appear after parents. Each pedigree occurs in a block
and is defined implicitly by the provided fathers and mothers.
Person i is a pedigree founder if and only if mother[i] = 0.
"""
function delta7_matrix(pedigree::Pedigree, person::Person,
  kinship::Matrix{Float64}, ped::Int)

  start = pedigree.start[ped]
  finish = pedigree.finish[ped]
  delta7 = zeros(pedigree.individuals[ped], pedigree.individuals[ped])
  q = start - 1
  for i = start:finish
    for j = start:i
      if j == i
        delta7[i - q, i - q] = 1.0
      elseif person.mother[i] > 0 && person.mother[j] > 0
        k = person.father[i]
        l = person.mother[i]
        m = person.father[j]
        n = person.mother[j]
        delta7[i - q, j - q] = kinship[k - q, m - q] * kinship[l - q, n - q] +
          kinship[k - q, n - q] * kinship[l - q, m - q]
        delta7[j - q, i - q] = delta7[i - q, j - q]
      end
    end
  end
  #
  # Extend the delta7 matrix to co-twins.
  #
  cotwin_extension!(delta7, pedigree, person, ped)
  return delta7
end # function delta7_matrix

"""
This function computes Jacquard's nine condensed identity
coefficients for each pedigree by simulation. Children
must appear after parents. A pedigree is defined implicitly
by giving fathers and mothers. Person i is a pedigree founder
if and only if mother[i] = 0. The more replicates requested,
the more precision attained. X-linked inheritance is possible.
"""
function jacquard_coefficients(pedigree::Pedigree, person::Person,
  ped::Int, replicates::Int, xlinked::Bool)

  start = pedigree.start[ped]
  finish = pedigree.finish[ped]
  ped_size = pedigree.individuals[ped]
  coefficient = zeros(9, ped_size, ped_size)
  x = zeros(ped_size, ped_size)
  source = zeros(Int, 2, ped_size)
  #
  # Loop over all replicates.
  #
  q = start - 1
  for rep = 1:replicates
    f = 0
    for j = start:finish
      if person.mother[j] == 0
        f = f + 1
        source[1, j - q] = f
        if xlinked && person.male[j]
        else
          f = f + 1
        end
        source[2, j - q] = f
      else
        if xlinked && person.male[j]
          u = rand(1:2, 1)
          source[1, j - q] = source[u[1], person.mother[j] - q]
          source[2, j - q] = source[1, j - q]
        else
          u = rand(1:2, 2)
          source[1, j - q] = source[u[1], person.father[j] - q]
          source[2, j - q] = source[u[2], person.mother[j] - q]
        end
      end
    end
    #
    # Loop over all pairs.
    #
    for j = start:finish
      for k = start:finish
        l = identity_state(source[:, j - q], source[:, k - q])
        coefficient[l, j - q, k - q] = coefficient[l, j - q, k - q] + 1.0
      end
    end
  end
  #
  # Extend the coefficients to co-twins.
  #
  for i = 1:9
    x = reshape(coefficient[i, :, :], ped_size, ped_size)
    cotwin_extension!(x, pedigree, person, ped)
    coefficient[i, :, :] = reshape(x, 1, ped_size, ped_size)
  end
  return coefficient / replicates
end # function jacquard_coefficients

"""
This function extends identity coefficients to co-twins.
"""
function cotwin_extension!(x::Matrix{Float64}, pedigree::Pedigree,
  person::Person, ped::Int)

  q = pedigree.start[ped] - 1
  for i = pedigree.finish[ped] + 1:pedigree.twin_finish[ped]
    m = person.primary_twin[i]
    for j = pedigree.start[ped]:pedigree.finish[ped]
      x[i - q, j - q] = x[m - q, j - q]
      x[j - q, i - q] = x[i - q, j - q]
    end
    for j = pedigree.finish[ped] + 1:pedigree.twin_finish[ped]
      n = person.primary_twin[j]
      x[i - q, j - q] = x[m - q, n - q]
    end
  end
end # function cotwin_extension!

"""
This function assigns one of the 9 condensed identity states
to two sourced genotypes. The integer n is the number of
unique genes among the 4 sampled genes.
"""
function identity_state(source1::Vector{Int}, source2::Vector{Int})

  n = length(union(IntSet(source1), IntSet(source2)))
  if n == 4
    return 9
  elseif n == 3
    if source1[1] == source1[2]
      return 4
    elseif source2[1] == source2[2]
      return 6
    else
      return 8
    end
  elseif n == 2
    if source1[1] == source1[2]
      if source2[1] == source2[2]
        return 2
      else
        return 3
      end
    else
      if source2[1] == source2[2]
        return 5
      else
        return 7
      end
    end
  else
    return 1
  end
end # function identity_state

"""Assembles names for printing."""

function make_plot_name(x::DataFrame)
  name = Array{String}(size(x, 1))
  for i in 1:length(name)
    name[i] = "Person1=" * x[i, 3] * ", " * "Person2=" * x[i, 4]
  end
  return name
end

"""Composes the comparison plot using PlotlyJS."""

function make_compare_plot(x::DataFrame, name::Vector{String})
  trace1 = scatter(;x=x[:theoretical_kinship],
    y=x[:empiric_kinship], mode="markers",
    name="empiric kinship", text=name)
#
  trace2 = scatter(;x=[1/2, 1/4, 1/8, 1/16, 1/32, 1/64, 1/128, 0.0],
    y=[1/2, 1/4, 1/8, 1/16, 1/32, 1/64, 1/128, 0.0],
    mode="markers", name="marker for midpoint")
#
  layout = Layout(;title="Compare empiric vs theoretical kinship",
    hovermode="closest",
    xaxis=attr(title="Theoretical kinship",
      showgrid=false, zeroline=false),
    yaxis=attr(title="Empiric Kinship", zeroline=false))
#
  data = [trace1, trace2]
  plot(data, layout)
end

"""Computes Fisher's atanh transform and standardizes it."""

function compute_fishers_z(x::DataFrame)
  theoretical_transformed = map(atanh, x[:theoretical_kinship])
  empiric_transformed = map(atanh, x[:empiric_kinship])
  difference = empiric_transformed - theoretical_transformed
  new_zscore = zscore(difference) #in StatsBase package
  return [x DataFrame(fishers_zscore = new_zscore)]
end

"""Plot the transformed data."""

function plot_fisher_z(x::DataFrame, name::Vector{String})
    trace1 = histogram(x=x[:fishers_zscore], text=name)
    data = [trace1]
#
    layout = Layout(barmode="overlay",
        title="Z-score plot for Fisher's statistic",
        xaxis=attr(title="Standard deviations"),
        yaxis=attr(title="count"))
#
    plot(data, layout)
end
##

end # module MendelKinship

max_stars_repo_name: chris-german/Tutorials

max_issues_repo_name: chris-german/Tutorials

max_forks_repo_name: chris-german/Tutorials


text:
get_pre(::Synapse, sys::ODESystem) = sys.V
get_post(::Synapse, sys::ODESystem) = sys.V

my_pre(::Synapse, syn_sys::AbstractTimeDependentSystem) = syn_sys.Vpre
my_post(::Synapse, syn_sys::AbstractTimeDependentSystem) = syn_sys.Vpost
post_connector(::Synapse) = :Isyn

mutable struct Chol{T<:AbstractFloat} <: Synapse
    ḡChol::T
    s::T
    Eₛ::T # mV
    k₋::T # ms
    Vth::T # mV
    δ::T   # mV
end

Chol(x) = Chol(x, 0.0, -80.0, 0.01, -35.0, 5.0)
s̄(syn::Chol, Vpre) = 1.0 / (1.0 + exp((syn.Vth - Vpre) / syn.δ))
τs(syn::Chol, Vpre) = (1.0 - s̄(syn, Vpre)) / syn.k₋
syn_current(::Chol, sys::ODESystem) = sys.IChol

mutable struct Glut{T<:AbstractFloat} <: Synapse
    ḡGlut::T
    s::T
    Eₛ::T
    k₋::T
    Vth::T
    δ::T
end

Glut(x) = Glut(x, 0.0, -70.0, 0.025, -35.0, 5.0)
s̄(syn::Glut, Vpre) = 1.0 / (1.0 + exp((syn.Vth - Vpre) / syn.δ))
τs(syn::Glut, Vpre) = (1.0 - s̄(syn, Vpre)) / syn.k₋
syn_current(::Glut, sys::ODESystem) = sys.IGlut

function channel_dynamics(ch::Chol, Vpre, Vpost)
    states = @variables s(t) IChol(t)
    parameters = @parameters ḡChol
    eqs = [D(s) ~ (1 / τs(ch, Vpre)) * (s̄(ch, Vpre) - s),
        IChol ~ -ḡChol * s * (Vpost - ch.Eₛ)]
    current = [eqs[2]]
    defaultmap = [s => ch.s, ḡChol => ch.ḡChol]
    return eqs, states, parameters, current, defaultmap
end

function channel_dynamics(ch::Glut, Vpre, Vpost)
    states = @variables s(t) IGlut(t)
    parameters = @parameters ḡGlut
    eqs = [D(s) ~ (1 / τs(ch, Vpre)) * (s̄(ch, Vpre) - s),
        IGlut ~ -ḡGlut * s * (Vpost - ch.Eₛ)]
    current = [eqs[2]]
    defaultmap = [s => ch.s, ḡGlut => ch.ḡGlut]
    return eqs, states, parameters, current, defaultmap
end

function (syn::Synapse)(; name=get_name(syn))
    @variables Vpre(t) Vpost(t)
    eqs, states, parameters, current, defaultmap = channel_dynamics(syn, Vpre, Vpost)
    sys = ODESystem(eqs, t, [Vpre, Vpost, states...], [parameters...]; observed=current,
        name=name, defaults=defaultmap)
    return sys
end

struct directed_synapse{S<:Synapse,N} <: Synapse
    pre_n::N
    post_n::N
    syn::S
end


text:
function value(obj::UnivariateMatrixFunction, arg0::jdouble)
    return jcall(obj, "value", Vector{Vector{jdouble}}, (jdouble,), arg0)
end


text:
#=
Given a list of integers, write a function that returns the largest sum of non-adjacent numbers. Numbers can be 0 or negative.

For example, [2, 4, 6, 2, 5] should return 13, since we pick 2, 6, and 5. [5, 1, 1, 5] should return 10, since we pick 5 and 5.

Follow-up: Can you do this in O(N) time and constant space?
=#

function sum_notadjacent(X::Array{Int, 1})
    sum = zeros(Int8, length(X))
    sum[1] = maximum([X[1], 0])
    if length(X) >= 2
        sum[2] = maximum([sum[1], X[2]])
    end
    if length(X) >= 3
        for index in 3:length(X)
            sum[index] = maximum([sum[index-2] + X[index], sum[index-1]])
        end
    end
    return sum[end]
end

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
#=========================================================================
This is a reimplementation of LinearAlgebra: rmul!, lmul!, ldiv!, rdiv!, where both arguments are triangular.
Note:
- Our triangular matrices can be any of the 4 subtypes of LinearAlgbra.AbstractTriangular, i.e. any of:
        {LowerTriangular, UpperTriangular, UnitLowerTriangular, UnitUpperTriangular}
- Adjoints and Transposes of triangular matrices are also treated as triangular.

This code addresses the issues raised here:
https://discourse.julialang.org/t/issues-with-inplace-updates-to-triangular-matrices-using-lmul-rmul-ldiv-and-rdiv/24887

This code uses nested loops to declare (using @eval) a multitude of very specific methods, thus shadowing
existing implementations in triangular.jl and bidiag.jl, which are problematic in V1.1.1.

Some properties of the new methods are:

- They don't call down into BLAS. All methods use three nested loops to iterate over martrix elements, using
  getindex! and setindex!.

- They are probably faster than pre-existing methods, for small matrices.

- For large matrices, the new methods can be slower or faster than pre-existing implementations, depending on whether
  the type combinations would have landed up (quite unpredictably) in triangular.jl or bidiag.jl.

- The new methods aim to be exhaustive for all combinations of arguments that would allow inplace triangle
  updates, i.e. where:
  -- The uplo of both arguments agree, and
  -- the element types agree.
  Note:
  -- The uplo of e.g. Transpose(LowerTriangular) and UnitUpperTriangular would agree.
  -- Transposes and Adjoints of triangular matrices can also serve as to-be-updated destinations.

- No clobbering of triangle parents outside of the intended triangle.

- For L,R triangular and fun! ∈ {rmul!, lmul!, ldiv!, rdiv!}:
      applicable(fun!,L,R) && fun!(L,R)
  should not crash, except when
  -- sizes disagree, or
  -- when divisors are singular, or
  -- when dividing matrices of integers (this should ideally be changed so that applicable returns false).

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

import LinearAlgebra: rmul!, lmul!, ldiv!, rdiv!

UpperA{T} = Adjoint{T,<:LowerTriangular{T}}
UpperT{T} = Transpose{T,<:LowerTriangular{T}}
upperset = (UpperTriangular,UpperA,UpperT)

UnitUpperA{T} = Adjoint{T,<:UnitLowerTriangular{T}}
UnitUpperT{T} = Transpose{T,<:UnitLowerTriangular{T}}
unit_upperset = (UnitUpperTriangular,UnitUpperT,UnitUpperA)

all_upperset = (upperset..., unit_upperset...)

for (DType,FType) in Iterators.product(all_upperset,all_upperset)
    if DType ∈ upperset # all cases where result is non-unit uppertriangular

        @eval begin

            # upper
            function rmul!(D::$DType{T}, F::$FType{T}) where {T}
                m,n = size(D,1), size(F,1)
                n == m || error("$m≠$n")
                @inbounds for j = n:-1:1
                    for i = j:-1:1
                        s = D[i,i]*F[i,j]
                        for k = i+1:j
                            s += D[i,k]*F[k,j]
                        end
                        D[i,j] = s
                    end
                end
                return D
            end

            # upper
            function lmul!(F::$FType{T}, D::$DType{T}) where {T}
                m,n = size(D,1), size(F,1)
                n == m || error("$m≠$n")
                @inbounds for j = 1:n
                    for i = 1:j
                        s = F[i,i]*D[i,j]
                        for k = i+1:j
                            s += F[i,k]*D[k,j]
                        end
                        D[i,j] = s
                    end
                end
                return D
            end

            # upper
            function ldiv!(F::$FType{T}, D::$DType{T}) where {T}
                m,n = size(D,1), size(F,1)
                n == m || error("$m≠$n")
                @inbounds for j = n:-1:1
                    (den = F[j,j]) == 0 && error("F is singular")
                    D[j,j] /= den
                    for i = j-1:-1:1
                        s = F[i,j] * D[j,j]
                        for k = i+1:j-1
                            s += F[i,k] * D[k,j]
                        end
                        (den = F[i,i]) == 0 && error("F is singular")
                        D[i,j] = (D[i,j] - s) / den
                    end
                end
                return D
            end

            # upper
            function rdiv!(D::$DType{T}, F::$FType{T}) where {T}
                m,n = size(D,1), size(F,1)
                n == m || error("$m≠$n")
                @inbounds for i = 1:n
                    (den = F[i,i]) == 0 && error("F is singular")
                    D[i,i] /= den
                    for j = i+1:n
                        s = D[i,i] * F[i,j]
                        for k = i+1:j-1
                            s += D[i,k] * F[k,j]
                        end
                        (den = F[j,j]) == 0 && error("F is singular")
                        D[i,j] = (D[i,j] - s) / den
                    end
                end
                return D
            end

        end # eval

    elseif DType ∈ unit_upperset && FType ∈ unit_upperset # all cases where result is unit uppertriangular

        @eval begin

            # unit upper
            function rmul!(D::$DType{T}, F::$FType{T}) where {T}
                m,n = size(D,1), size(F,1)
                n == m || error("$m≠$n")
                @inbounds for j = n:-1:1
                    for i = j-1:-1:1
                        s = F[i,j]
                        for k = i+1:j
                            s += D[i,k]*F[k,j]
                        end
                        D[i,j] = s
                    end
                end
                return D
            end

            # unit upper
            function lmul!(F::$FType{T}, D::$DType{T}) where {T}
                m,n = size(D,1), size(F,1)
                n == m || error("$m≠$n")
                @inbounds for j = 1:n
                    for i = 1:j-1
                        s = D[i,j]
                        for k = i+1:j
                            s += F[i,k]*D[k,j]
                        end
                        D[i,j] = s
                    end
                end
                return D
            end

            # unit upper
            function ldiv!(F::$FType{T}, D::$DType{T}) where {T}
                m,n = size(D,1), size(F,1)
                n == m || error("$m≠$n")
                @inbounds for j = n:-1:1
                    for i = j-1:-1:1
                        s = F[i,j]
                        for k = i+1:j-1
                            s += F[i,k] * D[k,j]
                        end
                        D[i,j] -= s
                    end
                end
                return D
            end

            # unit upper
            function rdiv!(D::$DType{T}, F::$FType{T}) where {T}
                m,n = size(D,1), size(F,1)
                n == m || error("$m≠$n")
                @inbounds for i = 1:n
                    for j = i+1:n
                        s = F[i,j]
                        for k = i+1:j-1
                            s += D[i,k] * F[k,j]
                        end
                        D[i,j] -= s
                    end
                end
                return D
            end

        end  #eval

    end # if
end # for

# lower triangular

LowerA{T} = Adjoint{T,<:UpperTriangular{T}}
LowerT{T} = Transpose{T,<:UpperTriangular{T}}
lowerset = (LowerTriangular,LowerA,LowerT)

UnitLowerA{T} = Adjoint{T,<:UnitUpperTriangular{T}}
UnitLowerT{T} = Transpose{T,<:UnitUpperTriangular{T}}
unit_lowerset = (UnitLowerTriangular,UnitLowerT,UnitLowerA)

all_lowerset = (lowerset..., unit_lowerset...)

for (DType,FType) in Iterators.product(all_lowerset,all_lowerset)
    if DType ∈ lowerset # all cases where result is non-unit lowertriangular

        @eval begin

            # lower
            function rmul!(D::$DType{T}, F::$FType{T}) where {T}
                m,n = size(D,1), size(F,1)
                n == m || error("$m≠$n")
                @inbounds for j=1:n
                    for i=j:n
                        s = D[i,j]*F[j,j]
                        for k=j+1:i
                            s += D[i,k]*F[k,j]
                        end
                        D[i,j] = s
                    end
                end
                return D
            end

            # lower
            function lmul!(F::$FType{T}, D::$DType{T}) where {T}
                m,n = size(D,1), size(F,1)
                n == m || error("$m≠$n")
                @inbounds for j=n:-1:1
                    for i=n:-1:j
                        s = F[i,j]*D[j,j]
                        for k=j+1:i
                            s += F[i,k]*D[k,j]
                        end
                        D[i,j] = s
                    end
                end
                return D
            end

            # lower
            function ldiv!(F::$FType{T}, D::$DType{T}) where {T}
                m,n = size(D,1), size(F,1)
                n == m || error("$m≠$n")
                @inbounds for j = 1:n
                    (den = F[j,j]) == 0 && error("F is singular")
                    D[j,j] /= den
                    for i = j+1:n
                        s = F[i,j] * D[j,j]
                        for k = j+1:i-1
                            s += F[i,k] * D[k,j]
                        end
                        (den = F[i,i]) == 0 && error("F is singular")
                        D[i,j] = (D[i,j] - s) / den
                    end
                end
                return D
            end

            # lower
            function rdiv!(D::$DType{T}, F::$FType{T}) where {T}
                m,n = size(D,1), size(F,1)
                n == m || error("$m≠$n")
                @inbounds for i = n:-1:1
                    (den = F[i,i]) == 0 && error("F is singular")
                    D[i,i] /= den
                    for j = i-1:-1:1
                        s = D[i,i] * F[i,j]
                        for k = j+1:i-1
                            s += D[i,k] * F[k,j]
                        end
                        (den = F[j,j]) == 0 && error("F is singular")
                        D[i,j] = (D[i,j] - s) / den
                    end
                end
                return D
            end

        end #eval

    elseif DType ∈ unit_lowerset && FType ∈ unit_lowerset # all cases where result is unit lowertriangular

        @eval begin

            # unit lower
            function rmul!(D::$DType{T}, F::$FType{T}) where {T}
                m,n = size(D,1), size(F,1)
                n == m || error("$m≠$n")
                @inbounds for j=1:n
                    for i=j+1:n
                        s = D[i,j]
                        for k=j+1:i
                            s += D[i,k]*F[k,j]
                        end
                        D[i,j] = s
                    end
                end
                return D
            end

            # unit lower
            function lmul!(F::$FType{T}, D::$DType{T}) where {T}
                m,n = size(D,1), size(F,1)
                n == m || error("$m≠$n")
                @inbounds for j=n:-1:1
                    for i=n:-1:j+1
                        s = F[i,j]
                        for k=j+1:i
                            s += F[i,k]*D[k,j]
                        end
                        D[i,j] = s
                    end
                end
                return D
            end

            # unit lower
            function ldiv!(F::$FType{T}, D::$DType{T}) where {T}
                m,n = size(D,1), size(F,1)
                n == m || error("$m≠$n")
                @inbounds for j = 1:n
                    for i = j+1:n
                        s = F[i,j]
                        for k = j+1:i-1
                            s += F[i,k] * D[k,j]
                        end
                        D[i,j] -= s
                    end
                end
                return D
            end

            # unit lower
            function rdiv!(D::$DType{T}, F::$FType{T}) where {T}
                m,n = size(D,1), size(F,1)
                n == m || error("$m≠$n")
                @inbounds for i = n:-1:1
                    for j = i-1:-1:1
                        s = F[i,j]
                        for k = j+1:i-1
                            s += D[i,k] * F[k,j]
                        end
                        D[i,j] -= s
                    end
                end
                return D
            end

        end  # eval

    end # if
end # for


text:
module TWA
    export calculate_distribution,average,variance,
    mcs_for_survival_probability,survival_probability,Weyl,mcs_chain,
    mcs_for_distributions,monte_carlo_integrate,MonteCarloSystem,coherent_Wigner_HWxSU2,coherent_Wigner_SU2,coherent_Wigner_HW

    using Distributed
    using Distributions
    using ProgressMeter
    using DiffEqCallbacks
    using SymEngine
    using ..ClassicalSystems
    using DataStructures
    using ..PhaseSpaces
    """
    ```julia
    struct PhaseSpaceDistribution
    ```
    This object represents a probability distribution in the phase space.
    Currently, the only  implementation is through [`coherent_Wigner_SU2`](@ref),
     [`coherent_Wigner_HW`](@ref), and [`coherent_Wigner_HWxSU2`](@ref)
    """
    struct PhaseSpaceDistribution
        probability_density
        sample
        ħ::Real
    end
    """
    ```julia
    struct MonteCarloSystem
    ```
    This object may be passed to [`monte_carlo_integrate`](@ref TWA.monte_carlo_integrate).
    Use [`mcs_for_averaging`](@ref TWA.mcs_for_averaging),
    [`mcs_for_variance`](@ref TWA.mcs_for_variance), and
    [`mcs_for_survival_probability`](@ref TWA.mcs_for_survival_probability) to generate them,
    and [`mcs_chain`](@ref TWA.mcs_chain) to join them together.
    """
    struct MonteCarloSystem
        f_initial
        f_loop
        reducer
        f_final
        N::Integer
        ts::AbstractVector{<:Real}
        distribution::PhaseSpaceDistribution
        function MonteCarloSystem(f_initial,f_loop,reducer,f_final,N::Integer,ts::Union{Nothing,AbstractVector{<:Real}},distribution::PhaseSpaceDistribution)
            if ts===nothing
                ts=[0.0]
            end
            return new(f_initial,f_loop,reducer,f_final,N,ts,distribution)
        end
    end

    function scale_to_index(v,r::LinRange)
            if v<r.start || r.stop <v
                return NaN
            end
            if r.len==1
                return 1
            end
            return Int(round(
                    ((v-r.start)/(r.stop-r.start))*(r.len-1)
                    +1))

        end
    """
    ```julia
    function monte_carlo_integrate(system::ClassicalSystems.ClassicalSystem,
        mc_system::MonteCarloSystem;
        tolerate_errors::Bool = true,
        show_progress::Bool = true,
        maxNBatch::Real=Inf,kargs...)
    ```
    This function is the backbone of this module. It performs a Monte Carlo integration-type procedure.
    The argument `mc_system` determines a `distribution::`[`PhaseSpaceDistribution`](@ref), an integer `N`, and a list of times `ts`.
    This function calls an initializing function also determined by `mc_system`. Then, it samples
    `N` random points `x` from `distribution`, which are integrated using the Hamiltonian given by `system`.
    For each trajectory and time in `ts`, an operation is performed, which is determined again by `mc_system`.
    A final operation is then performed and the result is returned. This may seem a bit abstract,
    but it is a very flexible system. The applications are more concrete; for example, if you
    generate `mc_system` using [`mcs_for_averaging(... observable=f ...)`](@ref TWA.mcs_for_averaging),
    then the initial operation is to generate an array of zeroes the same length as `ts`. Then, for each initial condition
    `x`, the result of `f(x(t))` is added to each element of the array. Finally, the array is overall divided by `N` and then
    returned. This is exactly a Monte Carlo integration of the function `f(x)` over the classical evolution
    of the distribution.

    This function uses all the available workers, but make sure to import this module in all of them.

    # Arguments
    - `system` should be an instance of [`ClassicalSystems.ClassicalSystem`](@ref).
    - `mc_system` should be an instance of  [`MonteCarloSystem`](@ref TWA.MonteCarloSystem).
    # Keyword arguments
    - `tolerate_errors` indicates that some errors in the integration may be ignored.
      This is useful because sometimes a one-in-a-million numerical instability may arise, and
      you may want to ignore it. If more than `100` errors occur consecutively, then then the
      procedure is stopped. Defaults to `true`.
    - `show_progress` is a boolean that toggles the progress bar. (Default is `true`).
    - `maxNBatch` is the maximum number of batch-sizes sent to each worker. Defaults to `inf`.
    - `kargs` are redirected to [`ClassicalSystems.integrate`].
    """
    function monte_carlo_integrate(system::ClassicalSystems.ClassicalSystem,
        mc_system::MonteCarloSystem;
        tolerate_errors::Bool = true,
        show_progress::Bool = true,
        maxNBatch::Real=Inf,kargs...)

        f_inicial=mc_system.f_initial
        f_loop=mc_system.f_loop
        reductora= mc_system.reducer
        finalizadora= mc_system.f_final

        distribution_sampler = mc_system.distribution.sample
        ts=mc_system.ts
        N=mc_system.N
        progressbar= nothing
        trabajadores=length(workers())
        if trabajadores==0
            trabajadores=1
        end
        Neff=min(Int(ceil(N/trabajadores)),maxNBatch)
        workerload=[]
        rest=N
        while rest>0
            load=Neff
            rest=rest-load
            extra=-min(0,rest)
            rest=max(rest,0)
            push!(workerload,load-extra)
        end
        channel = RemoteChannel(()->Channel{Int8}(N), 1)
        res=0

        workersready=0
        missedprogressloop=trabajadores
        @sync begin
            # this task prints the progress bar
            @async while true
                v= take!(channel)
                if v==0
                    if progressbar!=nothing
                        finish!(progressbar)
                    end
                    break
                elseif v==1
                    if workersready>=trabajadores
                        if progressbar==nothing && show_progress
                            progressbar=Progress(Int(N-missedprogressloop))
                        end
                        if show_progress
                            next!(progressbar)
                        end
                    else
                        missedprogressloop+=1
                    end

                else
                    workersready+=1
                end
            end

            # this task does the computation
            @async begin
               try
                   res=@distributed reductora for load in workerload
                        val=f_inicial()
                        ti=1
                        ti_error=1
                        nerrores=0
                        call_loop= (u,t,_)-> begin
                            if ti>=ti_error
                                f_loop(val,u,ti)
                            end
                            ti+=1
                        end

                        j=1
                        notimeevolution= (ts==[0.0])
                        cb=nothing
                        if !notimeevolution
                             cb=FunctionCallingCallback(call_loop;funcat=ts,func_start = false)
                        end
                        while j <=load
                            u₀=distribution_sampler()

                            ti=1
                            try
                                if notimeevolution
                                    call_loop(u₀,0.0,0)
                                else
                                    ClassicalSystems.integrate(system,u₀,last(ts);kargs...,save_everystep=false,callback=cb)
                                end
                                ti_error=1

                                put!(channel,  j==1 ? 2 : 1 )
                                j+=1
                                nerrores=0

                            catch e
                                if tolerate_errors==false
                                    error(e)
                                end
                                if nerrores>100
                                    println("Too many errors")
                                    error(e)

                                end
                                nerrores+=1
                                ti_error=ti
                            end

                        end
                        val
                    end

                finally
                    put!(channel, 0) # this tells the printing task to finish

                end
            end
        end

        return finalizadora(res)
    end

    """
    ```julia
    function mcs_chain(Fs::AbstractVector{MonteCarloSystem})
    ```
    Generates a [`MonteCarloSystem`](@ref) by chaining together those in the array `Fs = [monteCarloSystem1, monteCarloSystem2, ...]`.
    The output of the system generated is an array that has the outputs of each system, in the same order as `Fs`.

    If `Fs` may be any iterable. If it is not an `AbstractVector`, `collect` will be called on it.
    """
    function mcs_chain(Fs::AbstractVector{MonteCarloSystem})
        ts= Fs[1].ts
        N= Fs[1].N
        distribution=Fs[1].distribution
        if any(N!=F.N for F in Fs)
            error("All systems must have the same N")
        end
        if any(ts!=F.ts for F in Fs)
            error("All systems must have the same ts")
        end
        if any(distribution!=F.distribution for F in Fs)
            error("All systems must have the same phase space distribution")
        end
        f_inicial=function()
            [F.f_initial() for F in Fs]
        end
        f_loop=function(valores,u,t)
            for i in 1:length(Fs)
                Fs[i].f_loop(valores[i],u,t)
            end
        end
        reductora=function(last,next)
             [Fs[i].reducer(last[i],next[i]) for i in 1:length(Fs)]

        end
        f_final=function(res)
            [Fs[i].f_final(res[i]) for i in 1:length(Fs)]
        end
        return MonteCarloSystem(f_inicial,f_loop,reductora,f_final,N,ts,distribution)
    end
    """
    ```julia
    function mcs_chain(mcs1::MonteCarloSystem, mcs2::MonteCarloSystem, ...)
    ```
    Functionally equivalent to `mcs_chain([mcs1, mcs2, ...])`.
    """
    function mcs_chain(mcss::Vararg{MonteCarloSystem,N}) where {N}
        mcs_chain(mcss)
    end
    mcs_chain(mcss) = mcs_chain(collect(mcss))

    function function_from_expression(system::ClassicalSystems.ClassicalSystem,expression)
        varnames=ClassicalSystems.varnames(system)
        if typeof(expression)!=SymEngine.Basic && !isempty(methods(expression))#es una funcion
            if any(length(m.sig.parameters)-1==length(varnames) for m in methods(expression))
                return expression
            elseif any(length(m.sig.parameters)-1==1 for m in methods(expression))
                f(x...)=expression(x)
                return f
            else
                error("You should pass a function that takes a n-vector or n arguments, where n is the dimension of the phase space.")
            end

        end

        if !(Symbol.(SymEngine.free_symbols(SymEngine.Basic(expression))) ⊆ varnames)
            error("The expression $(expression) has unknown variables. Write it in terms only of $(ClassicalSystems.varnames(system))")
        end
        if  typeof(expression)==SymEngine.Basic
            expression=SymEngine.toString(expression) #esto no es lo mejor, seria mejor usar expresiones y no strings pero algo pasa con las raices cuadradas...

        end
        return Base.eval(Main,Meta.parse("($(join(varnames,',')))->($expression)"))
    end
    function functions_from_expressions(system::ClassicalSystems.ClassicalSystem,expressions)
        funciones_varianza=[function_from_expression(system,s) for s in expressions]
        return funciones_varianza
    end
    """
    ```julia
    function mcs_for_averaging(
        system::ClassicalSystems.ClassicalSystem;
        observable,
        ts::AbstractVector{<:Real}=[0.0],
        N::Integer,
        distribution::PhaseSpaceDistribution)
    ```
    Generates a [`MonteCarloSystem`](@ref) that computes
    ```math
        \\frac{1}{N} \\sum_{i=1}^N \\text{observable}(\\mathbf{x}_i(t))
    ```
    for each ``t`` in `ts`, where ``\\mathbf{x}_i`` will be sampled from `distribution`. The system produces an array
    the same size as `ts`, containing the result for each time.

    # Arguments
    - `system` should be an instance of [`ClassicalSystems.ClassicalSystem`](@ref),
    # Keyword arguments
    - `observable` can a function in the form `f(x::Vector)` or `f(x1,x2 ,..., x_n)` with
      n the dimension of the phase space determined by `system`. `observable` can also be
      an expression determining the operations between the `varnames` determined by system. For example,
      if `system` were a instance of [`ClassicalDicke.DickeSystem`](@ref Main.ClassicalDicke.DickeSystem), then `observable`
      could be `:(q+p^2 +Q)`, `f(Q, q, P, p) = q + p^2 - Q` or `f(x) = x[2] + x[4]^2 - x[1]`, which
      are all equivalent.
      See also the submodule [`Weyl`](@ref DickeModel.TWA.Weyl), which produces expressions corresponding to
      the [Weyl symbols](https://en.wikipedia.org/wiki/Wigner%E2%80%93Weyl_transform#The_inverse_map) of quantum observables.
      Note: `observable` can also be an array of observables, in which case an array is returned for each time.
    - `ts` should be a sorted array of times.
    - `N` is the number of points to sample. The bigger the more accurate.
    - `distribution` should be an instance of [`PhaseSpaceDistribution`](@ref)
    """
    function mcs_for_averaging(
        system::ClassicalSystems.ClassicalSystem;
        observable,
        ts::AbstractVector{<:Real}=[0.0],
        N::Integer,
        distribution::PhaseSpaceDistribution)
        onevar=false

        if !isa(observable,Array)
            observable=[observable]
            onevar=true
        end

        lengthexpr=length(observable)
        f_inicial= function()
            [zeros(lengthexpr) for i in 1:length(ts)]
        end
        varnames=ClassicalSystems.varnames(system)
        funciones_a_promediar=functions_from_expressions(system,observable)
        f_loop=function(valores,u,i)

            for fi in 1:lengthexpr
                valores[i][fi]+=funciones_a_promediar[fi](u...)
            end
        end
        f_final=function(res)
            res=res/N
            if onevar
                res=[re[1] for re in res]
            end
            return res
        end

        return MonteCarloSystem(f_inicial,f_loop,+,f_final,N,ts,distribution)
    end
    """
    ```julia
    function mcs_for_distributions(system::ClassicalSystems.ClassicalSystem;
        x,
        xs::Union{AbstractRange{<:Real},Nothing}=nothing,
        y=nothing,
        ts::Union{AbstractRange{<:Real},Nothing}=nothing,
        animate::Bool=false,
        ys::Union{AbstractRange{<:Real},Nothing}=nothing,
        N::Integer,
        distribution::PhaseSpaceDistribution)
    ```
    Generates a [`MonteCarloSystem`](@ref) that produces a multidimensional histogram
    allowing to visualize the expected value of observables under [`PhaseSpaceDistribution`](@ref) `distribution`.

    Several behaviors can be produced.
    - If `x` is `:t`, a matrix of dimensions `length(ts)` ``\\times`` `length(ys)`
      is produced. The coordinate `(tᵢ,yⱼ)` gives the probability of finding the observable
      `y` between `ys[j]` and `ys[j+1]` at the time `ts[i]`.
    - If `y` is `:t`, the same as above, changing `x` by `y`.
    - If neither `x` nor `y` are `:t`, and `animate = false`, then a matrix of dimensions `length(xs)` ``\\times`` `length(ys)`
      is produced.  The coordinate `(xᵢ,yⱼ)` gives the probability of finding the observable
      `x` between `xs[i]` and `xs[i+1]` and the observable
      `y` between `ys[j]` and `ys[j+1]`, averaging over all the times in `ts`
    - If neither `x` nor `y` are `t`, and `animate = true` then an array of matrices of dimensions `length(xs)` ``\\times`` `length(ys)`
      is produced, the same size as `ts`. For the `n`th matrix in this array, the coordinate `(xᵢ,yⱼ)` gives the probability of finding the observable
      `x` between `xs[i]` and `xs[i+1]` and the observable
      `y` between `ys[j]` and `ys[j+1]` at time `ts[n]`.

    # Arguments
    - `system` should be an instance of [`ClassicalSystems.ClassicalSystem`](@ref),
    # Keyword arguments
    - `x` and `y` can be `:t`, or an observable as described in the arguments of [`mcs_for_averaging`](@ref).
    - `xs` and `ys` should be given if  `x` and  `y` are not  `:t`, in wich case they should be
      range objects (e.g. `0:0.1:1`) containing the bins for the histogram.
    - `ts` should be a range object (e.g. `0:0.1:1`) of times.
    - `animate` should be a `Bool`, which determines the behavior if neither `x` nor `y` are `t`. Defaults to `true`
    - `N` is the number of points to sample. The bigger the more accurate.
    - `distribution` should be an instance of [`PhaseSpaceDistribution`](@ref)

    Note: If `ts` and `y` are both not passed, then `y` is set to `:t` and `ts` to `0:0`.
    """
    function mcs_for_distributions(system::ClassicalSystems.ClassicalSystem;x,xs::Union{AbstractRange{<:Real},Nothing}=nothing,
        y=nothing,ts::Union{AbstractRange{<:Real},Nothing}=nothing,animate::Bool=false,ys::Union{AbstractRange{<:Real},Nothing}=nothing,N::Integer,distribution::PhaseSpaceDistribution)
        onematrix=false
        #vars=[]
        if ts==nothing && y==nothing
            y=:t
            ts=0:0
        end
        if  ts==nothing
            error("You have to pass :ts")
        end

        if x==:t || y==:t || animate == false
            onematrix=true

            if x==:t
                xs=ts
            end
            if y==:t
                ys=ts
            end
        end
        xs=LinRange(xs)
        ys=LinRange(ys)
        f_inicial= function()
            if onematrix
                matrices=zeros(length(ys),length(xs))
            else
                matrices= [DefaultDict{Tuple{UInt16,UInt16},UInt32}(0) for i in ts]
            end
            return matrices
        end
        fx=nothing
        fy=nothing
        if x!=:t
            fx=functions_from_expressions(system,[x])[1]
        end
        if y!=:t
            fy=functions_from_expressions(system,[y])[1]
        end
        f_loop=nothing
        if onematrix==true
            if x==:t
                f_loop=function(matrices,u,i)
                    indu=scale_to_index(fy(u...),ys)
                    if !isnan(indu)
                        matrices[indu,i]+=1
                    end
                end
            elseif y==:t
                f_loop=function(matrices,u,i)
                    indu=scale_to_index(fx(u...),xs)
                    if !isnan(indu)
                        matrices[i,indu]+=1
                    end
                end
            else
                f_loop=function(matrices,u,i)
                    iy=scale_to_index(fy(u...),ys)
                    ix=scale_to_index(fx(u...),xs)
                    if !isnan(iy) && !isnan(ix)
                        matrices[iy,ix]+=1
                    end
                end
            end

        else
            f_loop=function(matrices,u,i)
                iy=scale_to_index(fy(u...),ys)
                ix=scale_to_index(fx(u...),xs)

                if !isnan(iy) && !isnan(ix)

                    matrices[i][iy,ix]+=1
                end
            end

        end
        if onematrix!=true

            datachannel = RemoteChannel(()->Channel{Any}(Inf), 1)

            resdat=[zeros(length(ys),length(xs)) for i in 1:length(ts)]
            flush_task=@async while true
                x= take!(datachannel)
                if x==nothing
                    close(datachannel)
                    break
                else
                    ra=rand(1:100)
                    for ri in length(resdat):-1:1
                        r=resdat[ri]
                        m=pop!(x)
                        for (k,v) in m
                            r[k...]+=v/N
                        end
                        yield()
                    end
                end
            end
       end

        function finalizacion(x)
            if onematrix==true
                return x/N
            end
            if length(workers())<=1 #if there is 1 worker, add does not get called
                add(x,nothing)
            end
            put!(datachannel,nothing)
            wait(flush_task)
            return resdat
        end
        function add(a,b)
            if onematrix==true
                a.+=b
                return a
            end

            for x in [a,b]
                if x===nothing
                    continue
                end

                put!(datachannel,x)
            end
            return nothing
        end
        return MonteCarloSystem(f_inicial,f_loop,add,finalizacion,N,ts,distribution)
    end
    """
    ```julia
    function calculate_distribution(system::ClassicalSystems.ClassicalSystem,
        [[[same kargs as mcs_for_distributions]]],
        kargs...)
    ```
    Calls [`mcs_for_distributions`](@ref) and then [`monte_carlo_integrate`](@ref) on the resulting `MonteCarloSystem`. Extra `kargs` are sent to the latter.
    """
    function calculate_distribution(system::ClassicalSystems.ClassicalSystem;x,N::Integer,
        animate::Bool=false,xs::Union{AbstractRange{<:Real},Nothing}=nothing,y=nothing,
        ts::Union{AbstractRange{<:Real},Nothing}=nothing,ys::Union{AbstractRange{<:Real},Nothing}=nothing,distribution::PhaseSpaceDistribution,kargs...)
        return monte_carlo_integrate(system,mcs_for_distributions(system;x=x,y=y,ts=ts,xs=xs,ys=ys,N=N,distribution=distribution,animate=animate);kargs...)
    end
    """
    ```julia
    function average(system::ClassicalSystems.ClassicalSystem,
        [[[same kargs as mcs_for_averaging]]],
        kargs...)
    ```
    Calls [`mcs_for_averaging`](@ref) and then [`monte_carlo_integrate`](@ref) on the resulting `MonteCarloSystem`. Extra `kargs` are sent to the latter.
    """
    function average(system::ClassicalSystems.ClassicalSystem;observable,N::Integer,ts::AbstractArray{<:Real}=[0.0],distribution::PhaseSpaceDistribution,kargs...)
        return monte_carlo_integrate(system,
            mcs_for_averaging(system;observable=observable,ts=ts,N=N,distribution=distribution);kargs...)
    end
    """
    ```julia
    function mcs_for_variance(
        system::ClassicalSystems.ClassicalSystem;
        observable,
        ts::AbstractVector{<:Real}=[0.0],
        N::Integer,
        distribution::PhaseSpaceDistribution,
        return_average::Bool = false)
    ```
    Generates a [`MonteCarloSystem`](@ref) that computes the variance
    ```math
        \\frac{1}{N} \\sum_{i=1}^N \\text{observable}^2(\\mathbf{x}_i(t)) - \\left (\\frac{1}{N}\\sum_{i=1}^N \\text{observable}(\\mathbf{x}_i(t))\\right )^2
    ```
    for each ``t`` in `ts`, where ``\\mathbf{x}_i`` will be sampled from `distribution`. The system produces an array
    the same size as `ts`, containing the result for each time.

    # Arguments
      - See arguments and keyword arguments that [`mcs_for_averaging`](@ref).
      - If `return_average = false` (default), only produces the variance, else it
        will produce a tuple, where the first element is the variance and the second
        the average.
    """
    function mcs_for_variance(system::ClassicalSystems.ClassicalSystem;observable,N::Integer,ts::AbstractArray{<:Real}=[0.0],
        return_average::Bool=false,distribution::PhaseSpaceDistribution)
        onevar=false

        if !isa(observable,Array)
            observable=[observable]
            onevar=true
        end
        tam=length(observable)

        apromediar=[observable;[(x)->f(x...)^2 for f in TWA.function_from_expression.((system,),observable)]]
        mcsa=TWA.mcs_for_averaging(system,observable=apromediar,ts=ts,N=N,distribution=distribution)
        mf_final=function(res)
            r=mcsa.f_final(res)
            vars=[[re[i+tam] - re[i]^2 for i in 1:tam] for re in r]
            if onevar
                vars=[v[1] for v in vars]
            end
            if return_average
                proms=[[re[i] for i in 1:tam] for re in r]
                if onevar
                    proms=[p[1] for p in proms]
                end
                vars=vars,proms
            end

            return vars
        end
        return MonteCarloSystem(mcsa.f_initial,mcsa.f_loop,+,mf_final,N,ts,distribution)
    end
    """
    ```julia
    function variance(system::ClassicalSystems.ClassicalSystem,
        [[[same kargs as mcs_for_variance]]],
        kargs...)
    ```
    Calls [`mcs_for_variance`](@ref) and then [`monte_carlo_integrate`](@ref) on the resulting `MonteCarloSystem`. Extra `kargs` are sent to the latter.
    """
    function variance(system::ClassicalSystems.ClassicalSystem;observable,N::Integer,ts::AbstractArray{<:Real}=[0.0],distribution::PhaseSpaceDistribution,return_average=false,kargs...)
        return TWA.monte_carlo_integrate(system,
            mcs_for_variance(system;observable=observable,ts=ts,N=N,distribution=distribution,return_average = return_average);kargs...)
    end
    """
    ```julia
    function mcs_for_survival_probability(
        system::ClassicalSystems.ClassicalSystem;
        N::Integer,
        ts::AbstractArray{<:Real},
        distribution::PhaseSpaceDistribution)
    ```
    Generates a [`MonteCarloSystem`](@ref) that computes the survival probability
    through Eq. (C.7) of Ref. [Villasenor2020](@cite) (with ``M=`` `N`, ``w=`` `distribution`).

    # Arguments
    - `system` should be an instance of [`ClassicalSystems.ClassicalSystem`](@ref),
    # Keyword arguments
    - `N` is the number of points to sample. The bigger the more accurate.
    - `ts` should be a sorted array of times.
    - `distribution` should be an instance of [`PhaseSpaceDistribution`](@ref)
    """
    function mcs_for_survival_probability(system::ClassicalSystems.ClassicalSystem;
        N::Integer,
        ts::AbstractArray{<:Real},
        distribution::PhaseSpaceDistribution)
        f_inicial= function()
            zeros(length(ts))
        end
        f_loop=function(valores,u,i)
                valores[i]+=distribution.probability_density(u)
        end
        f_final=function(res)
            L=length(ClassicalSystems.varnames(system))/2
            return res*(2*π*distribution.ħ)^L/N
        end

        return MonteCarloSystem(f_inicial,f_loop,+,f_final,N,ts,distribution)
    end
    """
    ```julia
    function survival_probability(system::ClassicalSystems.ClassicalSystem,
        [[[same kargs as mcs_for_survival_probability]]],
        kargs...)
    ```
    Calls [`mcs_for_survival_probability`](@ref) and then [`monte_carlo_integrate`](@ref) on the resulting `MonteCarloSystem`. Extra `kargs` are sent to the latter.
    """
    function survival_probability(system::ClassicalSystems.ClassicalSystem;distribution::PhaseSpaceDistribution,N::Integer,ts::AbstractArray{<:Real},kargs...)
        return TWA.monte_carlo_integrate(system,
            mcs_for_survival_probability(system;ts=ts,N=N,distribution=distribution);integate_backwards=true,kargs...)
    end
    """
    ```julia
    function coherent_Wigner_HW(;q₀::Real,p₀::Real,j::Real=1,ħ::Real=1/j)
    ```
    Returns a [`PhaseSpaceDistribution`](@ref) corresponding to the two-dimensional Wigner
    function of a coherent state of the Heisenberg-Weyl algebra (i.e. a standard coherent state)
    centered at `q₀,p₀`. A value of `ħ` may be passed, or pass `j` to set `ħ = 1/j`. (See Eq. (B.1) of Ref. [Villasenor2020](@cite))
    """
    function coherent_Wigner_HW(;q₀::Real,p₀::Real,j::Real=1,ħ::Real=1/j)
        σ=sqrt(ħ/2)
        gaussiana_qp=Distributions.MvNormal([q₀,p₀], σ)
        function W((q,p))
                return Distributions.pdf(gaussiana_qp,[q,p])
        end
        function sample()
            q,p=Distributions.rand(gaussiana_qp)
            return [q,p]
        end
        return PhaseSpaceDistribution(W,sample,ħ)
    end
    """
    ```julia
    function coherent_Wigner_SU2(;Q₀::Real,P₀::Real,j::Real=1,ħ::Real=1/j)
    ```
    Similar to [`coherent_Wigner_HW`](@ref), but for a coherent state for the SU(2)
    algebra (a Bloch coherent state).  The approximation given by Eq. (B.4) of Ref. [Villasenor2020](@cite) is used.
    """
    function coherent_Wigner_SU2(;Q₀::Real,P₀::Real,j::Real=1,ħ::Real=1/j)
        σ=sqrt(ħ/2)
        gaussiana_Θ=Distributions.Normal(0, σ)
        Rayleigh_Θ=Distributions.Rayleigh(σ)
        R(θ,φ,θ0,φ0)= #rotation that sends north pole to θ0,φ0 applied to a θ,φ
            (atan(sqrt((cos(θ0)*cos(φ)*sin(θ)-cos(θ)*sin(θ0))^2+sin(θ)^2*sin(φ)^2),cos(θ)*cos(θ0)+cos(φ)*sin(θ)*sin(θ0)),
            mod(atan(-cos(φ0)*sin(θ)*sin(φ)-cos(θ0)*cos(φ)*sin(θ)*sin(φ0)+cos(θ)*sin(θ0)*sin(φ0),-cos(θ0)*cos(φ)*cos(φ0)*sin(θ)+cos(θ)*cos(φ0)*sin(θ0)+sin(θ)*sin(φ)*sin(φ0)),2*pi))

        centroθ=PhaseSpaces.θ_of_QP(Q₀,P₀)
        centroφ=PhaseSpaces.φ_of_QP(Q₀,P₀)
        function Wθ(θ,φ)
            try
                Θ=acos(cos(θ)*cos(centroθ)+ sin(θ)*sin(centroθ)*cos(φ - centroφ))
                return Distributions.pdf(gaussiana_Θ,Θ)*sqrt(2*π)*σ *(1/(2*π*σ^2))
            catch
                return 0.0
            end
        end
        function W((Q,P))
            try

                 θ=PhaseSpaces.θ_of_QP(Q,P)
                 φ=PhaseSpaces.φ_of_QP(Q,P)
                 return Wθ(θ,φ)
            catch
                return 0.0
            end
        end

        function sample()
            Θ=Distributions.rand(Rayleigh_Θ)
            θ,φ = R(Θ,rand()*2*pi,centroθ,centroφ)
            return [PhaseSpaces.Q_of_θφ(θ,φ),PhaseSpaces.P_of_θφ(θ,φ)]
        end
        return PhaseSpaceDistribution(W,sample,ħ)
    end
    """
    ```julia
    function coherent_Wigner_HWxSU2(;Q₀::Real, q₀::Real,
                                    P₀::Real, p₀::Real,
                                    j::Real=1, ħ::Real=1/j)
    ```
    Produces the Wigner function corresponding to the tensor product of [`coherent_Wigner_HW`](@ref) and [`coherent_Wigner_SU2`](@ref).
    """
    function coherent_Wigner_HWxSU2(;Q₀::Real,q₀::Real,P₀::Real,p₀::Real,j::Real=1,ħ::Real=1/j)
        WSU2=coherent_Wigner_SU2(Q₀=Q₀,P₀=P₀,j=j,ħ=ħ)
        WHW=coherent_Wigner_HW(q₀=q₀,p₀=p₀,j=j,ħ=ħ)
        W((Q,q,P,p))=WSU2.probability_density((Q,P))*WHW.probability_density((q,p))
        function sample()
            Q,P=WSU2.sample()
            q,p=WHW.sample()
            return [Q,q,P,p]
        end
        return PhaseSpaceDistribution(W,sample,ħ)
    end
    """
    ```julia
    function coherent_Wigner_HWxSU2(u₀::AbstractVector{<:Real},j::Real=1,ħ::Real=1/j)
    ```
    Returns  `coherent_Wigner_HWxSU2(Q₀=u₀[1], q₀=u₀[2], P₀=u₀[3], p₀=u₀[4], j=j, ħ=ħ)`.
    """
    coherent_Wigner_HWxSU2(u₀::AbstractVector{<:Real};j::Real=1,ħ::Real=1/j)=coherent_Wigner_HWxSU2(Q₀=u₀[1],q₀=u₀[2],P₀=u₀[3],p₀=u₀[4],j=j,ħ=ħ)

    """
    ```julia
    function coherent_Wigner_SU2(u₀::AbstractVector{<:Real},j::Real=1,ħ::Real=1/j)
    ```
    Returns `coherent_Wigner_SU2(Q₀=u₀[1], P₀=u₀[2], j=j, ħ=ħ)`.
    """
    coherent_Wigner_SU2(u₀;j::Real=1,ħ::Real=1/j)=coherent_Wigner_SU2(Q₀=u₀[1],P₀=u₀[2],j=j,ħ=ħ)

    """
    ```julia
    function coherent_Wigner_HW(u₀::AbstractVector{<:Real},j::Real=1,ħ::Real=1/j)
    ```
    Returns `coherent_Wigner_HW(;q₀=u₀[1], p₀=u₀[2], j=j, ħ=ħ).
    """
    coherent_Wigner_HW(u₀;j::Real=1,ħ::Real=1/j)=coherent_Wigner_HW(q₀=u₀[1],p₀=u₀[2],j=j,ħ=ħ)

    module Weyl
        using SymEngine
        SymEngine.@vars q p
        """
        ```julia
        function n(j::Real)
        ```
        Returns the Weyl symbol of the number operator ``W(n̂)``, `(p²+  q² - ħ)/2`, where `ħ = 1/j`.
        """
        function n(j::Real)
            SymEngine.@vars q p
            return (j*(p^2 + q^2) - 1)/2 #Weyl(n̂)=(p²+q²- ħ)/2 y nuestras q,p están escaladas por sqrt(j)
        end
        """
        ```julia
        function n²(j::Real)
        ```
        Returns the Weyl symbol of the number operator squared ``W(n̂^2) = W(n̂)^2 - \\hbar^2/4``, where ``\\hbar = 1/j``. (Note that the extra term ``\\hbar^2/4``
        appears due to the non-commutativity of ``q̂`` and ``p̂``. See Ref. [Polkovnikov2013](@cite) for details on how to compute these expressions.
        """
        function n²(j::Real)
            return n(j)^2 - 1/4
        end

        #"The Moyal Representation for Spin" JOSEPH C. VÁRILLY AND JOSÉ M. GRACIA-BONDÍA Ann. of Phys. 190, 107-148 (1989)
        #https://doi.org/10.1016/0003-4916(89)90262-5
        macro SU2angles(ex)
            q= quote
                SymEngine.@vars Q P
                cosθ=(P^2+Q^2)/2-1
                sinθ=sqrt(1-cosθ^2)
                cosφ=Q/sqrt(P^2 + Q^2)
                sinφ=-P/sqrt(P^2 + Q^2)
                $ex
            end
            return esc(q)
        end
        """
        ```julia
        function Jz(j::Real)
        ```
        Returns the Weyl symbol of the operator ``\\hat{J}_z``. (See p. 114 of Ref. [Varilly1989](@cite))
        """
        function Jz(j::Real)
            return @SU2angles sqrt(j*(j+1))*cosθ
        end
        """
        ```julia
        function Jx(j::Real)
        ```
        Returns the Weyl symbol of the operator ``\\hat{J}_x``. (See p. 114 of Ref. [Varilly1989](@cite))
        """
        function Jx(j::Real)
            return @SU2angles sqrt(j*(j+1))*sinθ*cosφ
        end
        """
        ```julia
        function Jy(j::Real)
        ```
        Returns the Weyl symbol of the operator ``\\hat{J}_y``. (See p. 114 of Ref. [Varilly1989](@cite))
        """
        function Jy(j::Real)
            return @SU2angles sqrt(j*(j+1))*sinθ*sinφ
        end
        """
        ```julia
        function Jz²(j::Real)
        ```
        Returns the Weyl symbol of the operator ``\\hat{J}_z^2``. (See bottom of p. 128 of Ref. [Varilly1989](@cite))
        """
        function Jz²(j::Real)
            bj=sqrt(j*(j+1)*(2*j-1)*(2*j+3))
            return @SU2angles (j*(j+1))/3 + (bj/2)*(cosθ^2 - 1/3)
        end
        """
        ```julia
        function Jx²(j::Real)
        ```
        Returns the Weyl symbol of the operator ``\\hat{J}_x^2``. (See bottom of p. 128 of Ref. [Varilly1989](@cite))
        """
        function Jx²(j::Real)
            bj=sqrt(j*(j+1)*(2*j-1)*(2*j+3))
            return @SU2angles (j*(j+1))/3 + (bj/2)*((sinθ*cosφ)^2 - 1/3)
        end
        """
        ```julia
        function Jy²(j::Real)
        ```
        Returns the Weyl symbol of the operator ``\\hat{J}_y^2``. (See bottom of p. 128 of Ref. [Varilly1989](@cite))
        """
        function Jy²(j::Real)
            bj=sqrt(j*(j+1)*(2*j-1)*(2*j+3))
            return @SU2angles (j*(j+1))/3 + (bj/2)*((sinθ*sinφ)^2 - 1/3)
        end
    end
end


text:
using Pkg
Pkg.activate(@__DIR__)
using CSV
using DataFrames
using StatsBase
using RCall
using SigmaRidgeRegression
using LaTeXStrings
using LinearAlgebra
using MLJ
using Random
using ColorSchemes
using Plots
using StatsPlots
using PGFPlotsX

msd_filepath = joinpath(@__DIR__, "dataset", "YearPredictionMSD.txt")
msd = DataFrame(CSV.File(msd_filepath; header=false))

function feature_map(X)
    mean_features = X[:, 1:12]
    var_features = X[:, 13:24]
    sd_features = sqrt.(var_features)
    cov_features = X[:, 25:90]
    cor_features = zeros(size(cov_features))
    cnt = 0
    for offset in 1:11
        for i in 1:(12 - offset)
            cnt = cnt + 1
            cor_features[:, cnt] =
                cov_features[:, cnt] ./
                sqrt.(var_features[:, i] .* var_features[:, i + offset])
        end
    end
    grp = GroupedFeatures([12, 12, 66, 66])
    return (MLJ.table([mean_features sd_features cov_features cor_features]), grp)
end

Y = Float64.(msd[:, 1])
X, groups = feature_map(Matrix(msd[:, 2:91]));

train_idx = 1:463_715
test_idx = (1:51_630) .+ 463_715

loo_single_ridge = LooRidgeRegressor(;
    ridge=SingleGroupRidgeRegressor(; groups=groups, center=true, scale=true)
)

loo_sigmaridge = LooSigmaRidgeRegressor(; groups=groups, center=true, scale=true)

loo_multi_ridge = LooRidgeRegressor(;
    ridge=MultiGroupRidgeRegressor(; groups=groups, center=true, scale=true),
    rng=MersenneTwister(1),
)

holdout_glasso = TunedSeagull(; groups=groups, center=true, scale=true)

ns_subsample = [500; 1000; 5000]
n_montecarlo = 20
Random.seed!(10)

mse_array = Array{Float64}(undef, length(ns_subsample), n_montecarlo, 4)
time_array = Array{Float64}(undef, length(ns_subsample), n_montecarlo, 4)
λs_array = Array{Any}(undef, length(ns_subsample), n_montecarlo, 4)
#an initial run to make sure everything is precompiled followed by n_montecarlo runs
for n_montecarlo in [1; n_montecarlo]
for (k, n_subsample) in enumerate(ns_subsample)
    for j in Base.OneTo(n_montecarlo)
        train_idx_subsample = sample(train_idx, n_subsample; replace=false)
        resampling_idx = [(train_idx_subsample, test_idx)]
        for (i, mach) in
            enumerate([loo_sigmaridge, loo_single_ridge, loo_multi_ridge, holdout_glasso])
            time_array[k, j, i] = @elapsed begin
                _mach = machine(mach, X, Y)
                _eval = evaluate!(_mach; resampling=resampling_idx, measure=l2)
            end
            mse_array[k, j, i] = _eval.measurement[1]
            λs_array[k, j, i] = deepcopy(_mach.report.best_λs)
        end
    end
end
end

pgfplotsx()

push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\usepackage{amsmath}")
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\usepackage{amssymb}")
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\usepackage{bm}")
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\newcommand{\blambda}{\bm{\lambda}}")
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\newcommand{\risk}[1]{\bm{R}(#1)}")
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\usepackage[bbgreekl]{mathbbol}")
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\newcommand{\sigmacv}{\bbsigma}")

_orange = RGB{Float64}(0.933027, 0.665164, 0.198652)
_orange2 = RGB{Float64}(0.85004, 0.540122, 0.136212)

method_names =
    [L"\sigmacv-\textrm{Ridge}" L"\textrm{Single Ridge}" L"\textrm{Multi Ridge}" L"\textrm{Group Lasso}"]

_thickness_scaling = 1.8
mse_plot = plot(
    dotplot(
        method_names,
        mse_array[1, :, :];
        title=L"n\;=\;%$(ns_subsample[1])",
        frame=:box,
        grid=nothing,
        yguide="Mean squared error",
        label=nothing,
        markerstrokecolor=_orange2,
        markerstrokewidth=0.5,
        alpha=0.7,
        thickness_scaling=_thickness_scaling,
        ylim=(89, 110),
        color=_orange,
    ),
    dotplot(
        method_names,
        mse_array[2, :, :];
        title=L"n\;=\;%$(ns_subsample[2])",
        frame=:box,
        grid=nothing,
        label=nothing,
        markerstrokecolor=_orange2,
        markerstrokewidth=0.5,
        ylim=(89, 110),
        alpha=0.7,
        thickness_scaling=_thickness_scaling,
        color=_orange,
    ),
    dotplot(
        method_names,
        mse_array[3, :, :];
        title=L"n\;=\;%$(ns_subsample[3])",
        frame=:box,
        grid=nothing,
        label=nothing,
        markerstrokecolor=_orange2,
        markerstrokewidth=0.5,
        thickness_scaling=_thickness_scaling,
        alpha=0.7,
        ylim=(89, 110),
        color=_orange,
    );
    size=(1650, 400),
    layout=(1, 3),
)

savefig(mse_plot, "one_million_songs_mse.pdf")

time_plot = plot(
    dotplot(
        method_names,
        time_array[1, :, :];
        title=L"n\;=\;%$(ns_subsample[1])",
        frame=:box,
        grid=nothing,
        yguide="Time (seconds)",
        label=nothing,
        markerstrokecolor=_orange2,
        alpha=0.7,
        markerstrokewidth=0.5,
        yscale=:log10,
        ylim = (0.1, 100),
        yticks = ([0.1, 1,10,100], ["0.1","1","10","100"]),
        #ylim=(-0.5, 58.0),
        thickness_scaling=_thickness_scaling,
        color=_orange,
    ),
    dotplot(
        method_names,
        time_array[2, :, :];
        title=L"n\;=\;%$(ns_subsample[2])",
        frame=:box,
        grid=nothing,
        label=nothing,
        alpha=0.7,
        markerstrokecolor=_orange2,
        markerstrokewidth=0.5,
        thickness_scaling=_thickness_scaling,
        yscale=:log10,
        ylim = (0.1, 100),
        yticks = ([0.1, 1,10,100], ["0.1","1","10","100"]),
        #ylim=(-0.5, 58.0),
        color=_orange,
    ),
    dotplot(
        method_names,
        time_array[3, :, :];
        title=L"n\;=\;%$(ns_subsample[3])",
        frame=:box,
        grid=nothing,
        label=nothing,
        alpha=0.7,
        markerstrokecolor=_orange2,
        markerstrokewidth=0.5,
        thickness_scaling=_thickness_scaling,
        yscale=:log10,
        ylim = (0.1, 100),
        yticks = ([0.1, 1,10,100], ["0.1","1","10","100"]),
        #ylim=(-0.5, 58.0),
        color=_orange,
    );
    size=(1650, 400),
    layout=(1, 3),
)

savefig(time_plot, "one_million_songs_time.pdf")

_trunc = 20

λs_mean = min.(getindex.(λs_array, 1), 20)
λs_std = min.(getindex.(λs_array, 2), 20)
λs_cov = min.(getindex.(λs_array, 3), 20)
λs_cor = min.(getindex.(λs_array, 4), 20)

λs_sigma_ridge = [λs_mean[1, :, 1] λs_std[1, :, 1] λs_cov[1, :, 1] λs_cor[1, :, 1]]

feature_names = ["mean" "std" "cov" "cor"]
λs_names = [L"\hat{\lambda}_{\textrm{%$n}}" for n in feature_names]

λ_plot_params = (;
    frame=:box,
    grid=nothing,
    label="",
    markerstrokecolor=:purple,
    markerstrokewidth=0.5,
    markercolor=RGB{Float64}(205 / 256, 153 / 256, 255 / 256),
    thickness_scaling=1.7,
    ylim=(-0.9, _trunc + 0.8),
)

method_names_short =
    [L"\sigmacv-\textrm{Ridge}" L"\textrm{Single}" L"\textrm{Multi}" L"\textrm{GLasso}"]

plot(
    dotplot(
        method_names_short,
        λs_mean[1, :, :];
        yguide=L"\min\{\widehat{\lambda},20\}",
        title=λs_names[1],
    ),
    dotplot(method_names_short, λs_std[1, :, :]; title=λs_names[2]),
    dotplot(method_names_short, λs_cov[1, :, :]; title=λs_names[3]),
    dotplot(method_names_short, λs_cor[1, :, :]; title=λs_names[4]);
    size=(1500, 270),
    layout=(1, 4),
    λ_plot_params...
)
savefig("one_million_song_lambdas_n500.pdf")

plot(
    dotplot(method_names_short, λs_mean[2, :, :]; yguide=L"\min\{\widehat{\lambda},20\}"),
    dotplot(method_names_short, λs_std[2, :, :]),
    dotplot(method_names_short, λs_cov[2, :, :]),
    dotplot(method_names_short, λs_cor[2, :, :]);
    size=(1500, 260),
    layout=(1, 4),
    λ_plot_params...
)
savefig("one_million_song_lambdas_n1000.pdf")

plot(
    dotplot(method_names_short, λs_mean[3, :, :]; yguide=L"\min\{\widehat{\lambda},20\}"),
    dotplot(method_names_short, λs_std[3, :, :]),
    dotplot(method_names_short, λs_cov[3, :, :]),
    dotplot(method_names_short, λs_cor[3, :, :]);
    size=(1500, 260),
    layout=(1, 4),
    λ_plot_params...,
)

savefig("one_million_song_lambdas_n5000.pdf")


text:
# Measure the kinetic (one-body) and the potential (two-body) energy
function measure_energy(
    lattice::Lattice, UUp_Array::Array{T1,3}, UDn_Array::Array{T2,3},
    n̄_up::Array{T3,1}, n̄_dn::Array{T4,1}
) where {T1<:FloatType, T2<:FloatType, T3<:FloatType, T4<:FloatType}

    EOneBody, ETwoBody = 0, 0

    for i in eachindex(view(lattice.KMatrix, 1 : lattice.Ns, 1 : lattice.Ns))
        if lattice.KMatrix[i] != 0
            Dij_up = sum(UUp_Array[:, i[1], i[2]] .* n̄_up)
            Dij_dn = sum(UDn_Array[:, i[1], i[2]] .* n̄_dn)
            EOneBody += -lattice.t * (Dij_up + Dij_dn)
        end
    end

    for i = 1 : lattice.Ns
        Umatrix = UUp_Array[:, i, i] * UDn_Array[:, i, i]'
        nmatrix = n̄_up * n̄_dn'
        Diiii = sum(Umatrix .* nmatrix)
        ETwoBody += lattice.U * Diiii
    end

    real(EOneBody), real(ETwoBody)
end

function mc_integration(lattice::Lattice, qmc::QMC)

    Ek_Array = Vector{Float64}()    # kinetic energy
    Ep_Array = Vector{Float64}()    # potential energy
    E_Array = Vector{Float64}()     # total energy
    sign_Array = Vector{Float64}()  # sign
    wrapping_count = 0

    # initialize a random field configuration
    σ_field = 2 * (rand(Float64, (lattice.Ns, qmc.L)) .< 0.5) .- 1
    # initial propagation
    Q0_Up, D0_Up, T0_Up, Q0_Dn, D0_Dn, T0_Dn =
                                    matrix_product_stablized(σ_field, lattice, qmc)

    eigvals_up = eigvals(D0_Up * T0_Up * Q0_Up, sortby = abs)
    eigvals_dn = eigvals(D0_Dn * T0_Dn * Q0_Dn, sortby = abs)
    # recursion
    ZUp = main_recursion(lattice.Ns, lattice.nup, eigvals_up, false)
    ZDn = main_recursion(lattice.Ns, lattice.ndn, eigvals_dn, false)
    # projection
    #ZUp = particle_projection(lattice.Ns, lattice.nup, eigvals_up, false)
    #ZDn = particle_projection(lattice.Ns, lattice.ndn, eigvals_dn, false)

    for n = 1 : qmc.sample_size

        # sweep over the field
        for l = 1 : qmc.L
            # recalculate the decomposition over a certain period
            # to ensure numerical stability
            if wrapping_count == 15
                σ_field_shifted = circshift(σ_field, (0, -l + 1))
                Q0_Up, D0_Up, T0_Up, Q0_Dn, D0_Dn, T0_Dn =
                                    matrix_product_stablized(σ_field_shifted, lattice, qmc)
                wrapping_count = 0
            end
            # store the decomposition temporarily
            Q0_Up_temp, D0_Up_temp, T0_Up_temp = Q0_Up, D0_Up, T0_Up
            Q0_Dn_temp, D0_Dn_temp, T0_Dn_temp = Q0_Dn, D0_Dn, T0_Dn
            # (use the first time slice as an example)
            # the following lines convert B_L * ... * B_1 = Q_0 * D_0 * T_0 to
            # B_1 * B_L * ... * B_2 = B_1 * [Q_0 * D_0 * T_0 * B_1^-1] = B_1 * Q_l * D_l * T_l
            # which preserves the eigenvalues and the trace.
            σ = σ_field[:, l]
            Bl_Up, Bl_Dn = singlestep_matrix(σ, lattice, qmc)

            QRCP_update!(Q0_Up, D0_Up, T0_Up, inv(Bl_Up), 'R')
            QRCP_update!(Q0_Dn, D0_Dn, T0_Dn, inv(Bl_Dn), 'R')

            # loop over sites
            for i = 1 : lattice.Ns

                # propose a new configuration
                σ[i] *= -1
                Bl_Up, Bl_Dn = singlestep_matrix(σ, lattice, qmc)

                Ql_Up, Dl_Up, Tl_Up = QRCP_update!(copy(Q0_Up), copy(D0_Up), copy(T0_Up),
                                                    Bl_Up, 'L')
                Ql_Dn, Dl_Dn, Tl_Dn = QRCP_update!(copy(Q0_Dn), copy(D0_Dn), copy(T0_Dn),
                                                    Bl_Dn, 'L')

                eigvals_up = eigvals(Dl_Up * Tl_Up * Ql_Up, sortby = abs)
                eigvals_dn = eigvals(Dl_Dn * Tl_Dn * Ql_Dn, sortby = abs)
                # recursion
                ZUpNew = main_recursion(lattice.Ns, lattice.nup, eigvals_up, false)
                ZDnNew = main_recursion(lattice.Ns, lattice.ndn, eigvals_dn, false)
                # projection
                #ZUpNew = particle_projection(lattice.Ns, lattice.nup, eigvals_up, false)
                #ZDnNew = particle_projection(lattice.Ns, lattice.ndn, eigvals_dn, false)
                # caculate the acceptence ratio
                R = abs((ZUpNew / ZUp) * (ZDnNew / ZDn))
                accept_ratio = R / (1 + R)

                if rand() < accept_ratio

                    # accept, update the field
                    σ_field[i, l] *= -1
                    ZUp, ZDn = ZUpNew, ZDnNew
                    # store the decomposition temporarily
                    Q0_Up_temp, D0_Up_temp, T0_Up_temp = Ql_Up, Dl_Up, Tl_Up
                    Q0_Dn_temp, D0_Dn_temp, T0_Dn_temp = Ql_Dn, Dl_Dn, Tl_Dn

                else
                    # roll back the change
                    σ[i] *= -1
                end

                # update Q_0, D_0, T_0 at the last site
                if i == lattice.Ns
                    Q0_Up, D0_Up, T0_Up = Q0_Up_temp, D0_Up_temp, T0_Up_temp
                    Q0_Dn, D0_Dn, T0_Dn = Q0_Dn_temp, D0_Dn_temp, T0_Dn_temp
                end

            end

            wrapping_count += 1

        end

        # measure observables
        # propagate using normal order
        MatProdUp_Q, MatProdUp_D, MatProdUp_T, MatProdDn_Q, MatProdDn_D, MatProdDn_T =
                                        matrix_product_stablized(σ_field, lattice, qmc)
        eigen_up = eigen(MatProdUp_D * MatProdUp_T * MatProdUp_Q, sortby = abs)
        eigen_dn = eigen(MatProdDn_D * MatProdDn_T * MatProdDn_Q, sortby = abs)
        # transformation matrix (A = P D P^-1)
        MatProdUp_P = MatProdUp_Q * eigen_up.vectors
        MatProdDn_P = MatProdDn_Q * eigen_dn.vectors
        # overlap coefficient
        UUp_Array, UDn_Array = overlap_coefficient(lattice.Ns, MatProdUp_P, MatProdDn_P)
        # recursion
        ZUp, n̄up = main_recursion(lattice.Ns, lattice.nup, eigen_up.values, true)
        ZDn, n̄dn = main_recursion(lattice.Ns, lattice.ndn, eigen_dn.values, true)
        # projection
        #ZUp, n̄up = particle_projection(lattice.Ns, lattice.nup, eigen_up.values, true)
        #ZDn, n̄dn = particle_projection(lattice.Ns, lattice.ndn, eigen_dn.values, true)
        # sign
        sgn = sign(ZUp) * sign(ZDn)
        # measure energy
        Ek, Ep = measure_energy(lattice, UUp_Array, UDn_Array, n̄up, n̄dn)
        # store data
        push!(sign_Array, sgn)
        push!(Ek_Array, Ek * sgn)
        push!(Ep_Array, Ep * sgn)
        push!(E_Array, (Ek + Ep) * sgn)

    end

    # block average
    EMean = mean(E_Array)
    EVar = varm(E_Array, EMean) / qmc.sample_size
    EkMean = mean(Ek_Array)
    EkVar = varm(Ek_Array, EkMean) / qmc.sample_size
    EpMean = mean(Ep_Array)
    EpVar = varm(Ep_Array, EpMean) / qmc.sample_size
    avgsgn = mean(sign_Array)

    EMean, EVar, EkMean, EkVar, EpMean, EpVar, avgsgn

end

max_stars_repo_name: TongSericus/CanonicalAFQMC

max_issues_repo_name: TongSericus/CanonicalAFQMC

max_forks_repo_name: TongSericus/CanonicalAFQMC


text:
using Blink
using Test
using Sockets

cleanup = !AtomShell.isinstalled()

cleanup && AtomShell.install()

@testset "basic functionality" begin
    # open window and wait for it to initialize
    w = Window(Blink.@d(:show => false), async=false);

    # make sure the window is really active
    @test @js(w, Math.log(10)) ≈ log(10)

    @test string(Blink.jsstring(:(Dict("a" => 1, :b => 10)))...) == "{\"a\":1,\"b\":10}"

    # check that <!DOCTYPE html> was declared
    @test  startswith(Blink.maintp.tokens[1].value, "<!DOCTYPE html>\n")
end

include("content/api.jl");
include("AtomShell/window.jl");

if Sys.iswindows()
    # Uninstalling AtomShell on Windows is currently broken:
    # https://github.com/JunoLab/Blink.jl/pull/143#issuecomment-414144008
    cleanup = false
end
cleanup && AtomShell.uninstall()


text:
using MultivariateStats
using Test
using LinearAlgebra
using StableRNGs

@testset "Regression" begin

    rng = StableRNG(34568)

    ## data

    m = 9
    n = 6
    n2 = 3

    X = randn(rng, m, n)
    A = randn(rng, n, n2)
    Xt = X'

    b = randn(rng, 1, n2)

    E = randn(rng, m, n2) * 0.1
    Y0 = X * A + E
    Y1 = X * A .+ b + E

    y0 = Y0[:,1]
    y1 = Y1[:,1]

    ## llsq

    A = llsq(X, Y0; dims=1, bias=false)
    A_r = copy(A)
    @test size(A) == (n, n2)
    @test X'Y0 ≈ X'X * A

    a = llsq(X, y0; dims=1, bias=false)
    @test size(a) == (n,)
    @test a ≈ A[:,1]

    A = llsq(Xt, Y0; dims=2, bias=false)
    @test size(A) == (n, n2)
    @test A ≈ A_r

    a = llsq(Xt, y0; dims=2, bias=false)
    @test size(a) == (n,)
    @test a ≈ A[:,1]

    Aa = llsq(X, Y1; dims=1, bias=true)
    Aa_r = copy(Aa)
    @test size(Aa) == (n+1, n2)
    A, b = Aa[1:end-1,:], Aa[end:end,:]
    @test X' * (Y1 .- b) ≈ X'X * A

    aa = llsq(X, y1; dims=1, bias=true)
    @test aa ≈ Aa[:,1]

    Aa = llsq(Xt, Y1; dims=2, bias=true)
    @test Aa ≈ Aa_r

    aa = llsq(Xt, y1; dims=2, bias=true)
    @test aa ≈ Aa[:,1]

    # test default dims=2 argument
    aa = llsq(X, y1)
    @test aa ≈ Aa[:,1]
    @test_throws DimensionMismatch llsq(X, y1; dims=2)
    @test_throws DimensionMismatch llsq(Xt, y1)

    ## ridge (with Real r)

    r = 0.1

    A = ridge(X, Y0, r; dims=1, bias=false)
    A_r = copy(A)
    @test size(A) == (n, n2)
    @test X'Y0 ≈ (X'X + r * I) * A

    a = ridge(X, y0, r; dims=1, bias=false)
    @test size(a) == (n,)
    @test a ≈ A[:,1]

    A = ridge(Xt, Y0, r; dims=2, bias=false)
    @test size(A) == (n, n2)
    @test A ≈ A_r

    a = ridge(Xt, y0, r; dims=2, bias=false)
    @test size(a) == (n,)
    @test a ≈ A[:,1]

    Aa = ridge(X, Y1, r; dims=1, bias=true)
    Aa_r = copy(Aa)
    @test size(Aa) == (n+1, n2)
    A, b = Aa[1:end-1,:], Aa[end:end,:]
    @test X' * (Y1 .- b) ≈ (X'X + r * I) * A

    aa = ridge(X, y1, r; dims=1, bias=true)
    @test aa ≈ Aa[:,1]

    Aa = ridge(Xt, Y1, r; dims=2, bias=true)
    @test Aa ≈ Aa_r

    aa = ridge(Xt, y1, r; dims=2, bias=true)
    @test aa ≈ Aa[:,1]

    # test default dims=2 argument
    aa = ridge(X, y1, r)
    @test aa ≈ Aa[:,1]
    @test_throws DimensionMismatch ridge(X, y1, r; dims=2)
    @test_throws DimensionMismatch ridge(Xt, y1, r)

    ## ridge (with diagonal r)

    r = 0.05 .+ 0.1 .* rand(rng, n)

    A = ridge(X, Y0, r; dims=1, bias=false)
    A_r = copy(A)
    @test size(A) == (n, n2)
    @test X'Y0 ≈ (X'X + diagm(0=>r)) * A

    a = ridge(X, y0, r; dims=1, bias=false)
    @test size(a) == (n,)
    @test a ≈ A[:,1]

    A = ridge(Xt, Y0, r; dims=2, bias=false)
    @test size(A) == (n, n2)
    @test A ≈ A_r

    a = ridge(Xt, y0, r; dims=2, bias=false)
    @test size(a) == (n,)
    @test a ≈ A[:,1]

    Aa = ridge(X, Y1, r; dims=1, bias=true)
    Aa_r = copy(Aa)
    @test size(Aa) == (n+1, n2)
    A, b = Aa[1:end-1,:], Aa[end:end,:]
    @test X' * (Y1 .- b) ≈ (X'X + diagm(0=>r)) * A

    aa = ridge(X, y1, r; dims=1, bias=true)
    @test aa ≈ Aa[:,1]

    Aa = ridge(Xt, Y1, r; dims=2, bias=true)
    @test Aa ≈ Aa_r

    aa = ridge(Xt, y1, r; dims=2, bias=true)
    @test aa ≈ Aa[:,1]

    ## ridge (with qudratic r matrix)

    Q = qr(randn(rng, n, n)).Q
    r = Q' * diagm(0=>r) * Q

    A = ridge(X, Y0, r; dims=1, bias=false)
    A_r = copy(A)
    @test size(A) == (n, n2)
    @test X'Y0 ≈ (X'X + r) * A

    a = ridge(X, y0, r; dims=1, bias=false)
    @test size(a) == (n,)
    @test a ≈ A[:,1]

    A = ridge(Xt, Y0, r; dims=2, bias=false)
    @test size(A) == (n, n2)
    @test A ≈ A_r

    a = ridge(Xt, y0, r; dims=2, bias=false)
    @test size(a) == (n,)
    @test a ≈ A[:,1]

    Aa = ridge(X, Y1, r; dims=1, bias=true)
    Aa_r = copy(Aa)
    @test size(Aa) == (n+1, n2)
    A, b = Aa[1:end-1,:], Aa[end:end,:]
    @test X' * (Y1 .- b) ≈ (X'X + r) * A

    aa = ridge(X, y1, r; dims=1, bias=true)
    @test aa ≈ Aa[:,1]

    Aa = ridge(Xt, Y1, r; dims=2, bias=true)
    @test Aa ≈ Aa_r

    aa = ridge(Xt, y1, r; dims=2, bias=true)
    @test aa ≈ Aa[:,1]

    ## isotonic
    xx = [3.3, 3.3, 3.3, 6, 7.5, 7.5]
    yy = [4, 5, 1, 6, 8, 7.0]
    a = isotonic(xx, yy)
    b = [[1,2,3],[4],[5,6]]
    @test a ≈ xx atol=0.1

    ## using various types
    @testset for (T, TR) in [(Int,Float64), (Float32, Float32)]
        X, y = rand(T, 10, 3), rand(T, 10)
        a = llsq(X, y)
        @test eltype(a) == TR
        a = ridge(X, y, one(T))
        @test eltype(a) == TR
    end

end


text:
# ######################
# A DataFrame way to get file

import Base: getindex, view

export getindex

Base.getindex(file::JDFFile, rows, col::Symbol) = begin
	loadjdf(file; cols = [col])[rows, 1]
end

Base.getindex(file::JDFFile, rows, cols::AbstractVector{Symbol}) = begin
	loadjdf(file; cols = cols)[rows, :]
end

Base.view(file::JDFFile, rows, cols) = getindex(file, rows, cols)

getindex(file::JDFFile, rows, cols) = loadjdf(file)[rows, cols]


text:
using KernelDensity
using Random
using QuadGK
Random.seed!(1)

Zs = StandardNormalSample.(randn(2000) .* sqrt(2))

kernel1 = Empirikos.FlatTopKernel(1.0)
kernel2 = Empirikos.FlatTopKernel(0.2)

@test quadgk(z->pdf(kernel1, z), -50.0, 50.0)[1] ≈ 1.0 atol = 0.005
@test quadgk(z->pdf(kernel2, z), -50.0, 50.0)[1] ≈ 1.0 atol = 0.005

kde1 = kde(response.(Zs), kernel2; boundary = (-12.0,12.0), npoints=2048*4)

idxs = [2000; 2048*2; 6000]
for idx in idxs
    x1 = kde1.x[idx]
    kde1_f = kde1.density[idx]
    kde1_manual = mean( pdf.(kernel2, response.(Zs) .- x1) )
    @test kde1_f ≈ kde1_manual atol = 0.0001
end

# check the full fit using InfinityNormDensityBand
# at α = 0.5 and α = 0.2
Zs_half_std = NormalSample.(response.(Zs), 0.5)

floc_02 = InfinityNormDensityBand(a_min=-3.0, a_max=3.0, α =0.2 )
floc_05 = InfinityNormDensityBand(a_min=-3.0, a_max=3.0, α =0.5 )

floc_fit_02 = fit(floc_02, Zs)
floc_fit_05 = fit(floc_05, Zs)

floc_fit_half_02 = fit(floc_02, Zs_half_std)

@test std(floc_fit_02.midpoints[1]) == 1.0
@test std(floc_fit_half_02.midpoints[1]) == 0.5

# check default bandwidth is what we expect

@test floc_fit_02.method.kernel.h == 1/sqrt(log(2000))
@test floc_fit_05.method.kernel.h == floc_fit_02.method.kernel.h
@test floc_fit_half_02.method.kernel.h == 0.5/sqrt(log(2000))

# check c_infty_boot are the same in the two cases
@test floc_fit_02.boot_samples == floc_fit_05.boot_samples
@test floc_fit_05.C∞ == median(floc_fit_05.boot_samples)
@test floc_fit_02.C∞  == quantile(floc_fit_02.boot_samples, 0.8)
@test floc_fit_02.C∞ > floc_fit_05.C∞

# double check KDE fits using explicit formula as above
@test floc_fit_02.estimated_density == floc_fit_05.estimated_density
@test floc_fit_02.midpoints == floc_fit_05.midpoints

@test length(floc_fit_02.midpoints) == length(floc_fit_02.estimated_density)
rand_idxs = sample(1:length(floc_fit_02.midpoints), 20)

_kernel = floc_fit_02.method.kernel

for idx in rand_idxs
    loc = response(floc_fit_02.midpoints[idx])
    dens = floc_fit_02.estimated_density[idx]
    dens_manual = mean( pdf.(_kernel, response.(Zs) .- loc) )
    @test dens ≈ dens_manual atol = 0.0001
end