text:
module OncoDiversity

using DataFrames
using CSV
using Query
using Statistics
using Dates
import DataFrames: rename!

export diversity, calcOverQ, findInflection, findInflectionLocal, reportInflection, calcIPdiversity, patientdiversity,
  rename!, requirecol, countframe, countframefull, summarize, filteredDataset,
  Dataset, CDR3SeqData,
  AbstractAnalysis, CDR3SeqAnalysis,
  DiversityScores

abstract type AbstractAnalysis end
struct CDR3SeqAnalysis <: AbstractAnalysis end
struct SCRNASeqAnalysis<: AbstractAnalysis end
struct FlowCytoAnalysis <: AbstractAnalysis end
struct CytokineAnalysis <: AbstractAnalysis end
struct SomaScanAnalysis <: AbstractAnalysis end

abstract type Dataset end

"""    CDR3SeqData{T,U}

store two data tables that make up a CDR3 Dataset,
the per (patient, sequence) count and the physiochemisty data for each sequence.
These tables can be joined together usind `DataFrame(ds::CDR3SeqData)`.

"""
struct CDR3SeqData{T,U} <: Dataset
  counts::T
  physiochem::U
end

"""    diversity(seq, q)

Computes the generalized diversity index for a given value of q and the associated frequencies of the dataset.

```math
  ^qD = \\left( \\sum_i^n (p_i^q) \\right)^{1/(1-q)}
```
where p_i are the frequency across types and q is the order of diversity and you are summing over all types present in the dataset.

"""
function diversity(seq, q)
  0 <= q || error("q=$(q) must be greater than zero")

  if(q == 1.0)
    return exp(-sum(seq .* log.(seq)))
  else
    return sum(seq.^q)^(1/(1-q))
  end

end

"""    calcOverQ(origDF::DataFrame,rangeQ::Array{Float64})

Computes the diversity score over a range of q values. Returns a DataFrame where each column is a different q value and each row is a different patient.

Input: DataFrame of frequencies per patient and range of q
Output: DataFrame of diversity scores across the range of q per patient

"""
function calcOverQ(origDF::DataFrame,rangeQ::Array{Float64})
  qDdf = DataFrame(patient=unique(origDF[:,:patient]))

  for q in rangeQ
      tempDF = patientdiversity(origDF,q);
      qDdf[:,Symbol(q)]=tempDF[:,:diversity]
  end

  return qDdf
end

"""    findInflection(q::Array, qD::Array)

Computes the inflection point of a sequence numerically.

"""
function findInflection(q::Array, qD::Array)
  approxLogDeriv = q[2:end].*diff(qD)./diff(q);
  return approxLogDeriv;
end

"""    reportInflection(q::Array, qD::Array)

Computes the inflection point and the magnitude of the slope at that point using findInflection and findInflectionLocal functions.

"""
function reportInflection(q::Array, qD::Array)
  diffVector = findInflection(q, qD)
  inflectPt = argmin(diffVector)
  slopeInflectPt = abs(diffVector[inflectPt])
  qInflecPt = q[inflectPt+1]
  return (qInflecPt, slopeInflectPt)
end

"""    calcIPdiversity(df::DataFrame, qrange::Array)

Computes diversity over a range of q values (qrange) and an input DataFrame of frequencies across types, then numerically determines the inflection point and the magnitude of the slope at that point, and returns a DataFrame with a column for patient name, a column for q at which the inflection point occurs, and a column for the magnitude of the slope at the inflection point.

Input: DataFrame and array of q values
Output: DataFrame of inflection point q and slope diversity per patient

"""
function calcIPdiversity(df::DataFrame, qrange::Array)

  # set up new dataframe to capture the inflection point (InfPt) and slope at the inflection point (InflPtS)
  newDF = DataFrame(
            patient=unique(df[:,:patient]),
            InflPt=zeros(size(unique(df[:,:patient]))[1]),
            InflPtS=zeros(size(unique(df[:,:patient]))[1]))

  # calculate qD over a supplied range of q's
  dfRangeQ = calcOverQ(df,qrange)

  # calculate the IP & slope and fill into the dataframe
  counter = 1;
  for r in eachrow(dfRangeQ)

      qDs = collect(r[2:end])

      (IP,IPslope) = reportInflection(qrange,qDs)
      newDF[counter,:InflPt] = IP
      newDF[counter,:InflPtS] = IPslope
      counter = counter+1;
  end

  return newDF
end

"""    patientdiversity(countframe::DataFrame,q::Real)

Computes a DataFrame of diversity scores across patients from a countframe DataFrame and a single value of q.

Input: DataFrame and a real number
Output: DataFrame

"""
function patientdiversity(countframe::DataFrame,q::Real)
  #@show totUniq = length(unique(df.CDR3))

  diversityframe = @from i in countframe begin
          @group i by i.patient into g
          @select {
              patient=first(g.patient),
              diversity=diversity(g.count ./ sum(g.count), q),
              numclust = length(g.count),
              checkrange = length(g.count) > diversity(g.count ./ sum(g.count), q)
          }
      @collect DataFrame
  end
  return diversityframe
end

"""    requirecol(df::DataFrame, colname::AbstractString)

Raises an assertion error if colname is not a column of DataFrame.

"""
function requirecol(df::DataFrame, colname::AbstractString)
  @assert colname in names(df) "Expected column named: $colname, names are $(names(df))"
end

rename!(df::DataFrame, s::String, t::Symbol) = rename!(df, Symbol(s)=>t)

"""    countframe(a::CDR3SeqAnalysis, typeDF::DataFrame)

Counts the number of occurances of CDR3 by patients, after checking that patient and CDR3 columns are present in the DataFrame, and returns as a DataFrame of the counts.

Input: CDR3SeqAnalysis struct and DataFrame of CDR3 recoveries
Output: DataFrame

"""
function countframe(a::CDR3SeqAnalysis, typeDF::DataFrame)
  ## checks that you have a column of individual identifies and CDR3 sequences
  requirecol(typeDF, "patient")
  requirecol(typeDF, "CDR3")

  ## check if there is a column called abundances -- for cases where you do not need to count individual reads

  ## try seeing if there is a column called abundance
  ## otherwise assumes that 1 row in the dataframe = 1 read
  try
    requirecol(typeDF, "abundance")
    ## if requirecol not satisfied, will jump to catch

    countDF = @from i in typeDF begin
      @where i.abundance > 0
      @select {
        patient=i.patient,
        CDR3=i.CDR3,
        count=i.abundance
      }
      @collect DataFrame
    end;

  catch
    countDF = @from i in typeDF begin
      @group i by i.patient, i.CDR3 into g
      @select {
          patient=first(g.patient),
          CDR3=first(g.CDR3),
          count=length(g.CDR3)
      }
      @collect DataFrame
    end;
  end

  #return countDF

end

"""    countframefull(a::CDR3SeqAnalysis, df::DataFrame)

Merges df, a DataFrame of counts, with all the physiochemical information associated with the CDR3 sequences.

Input: CDR3SeqAnalysis struc and DataFrame
Output: DataFrame

"""
function countframefull(a::CDR3SeqAnalysis, df::DataFrame)
  rename!(df, "Filename", :patient)
  return innerjoin(countframe(CDR3SeqAnalysis(), df), df, on=[:patient, :CDR3])
end

function CDR3SeqData(df::DataFrame)
  requirecol(df, "CDR3")
  try
    requirecol(df, "patient")
  catch AssertionError
    requirecol(df, "Filename")
    rename!(df, "Filename", :patient)
  end
  return CDR3SeqData(
    countframe(CDR3SeqAnalysis(), df),
    select(df, Not(:patient)) |> @groupby(_.CDR3) |> x-> map(first, x) |> DataFrame
  )
end

DataFrame(data::CDR3SeqData) = innerjoin(data.counts, data.physiochem, on=[:CDR3])

"""    summarize(df::CDR3SeqData)

Computes a summary of the CDR3 sequences recovered in the dataset per patient including the total number of recoveries, minimum number on recoveries per sequence, maximum number of recoveries per sequence, average number of recoveries per sequence, median number of recoveries per sequence, and number of unique CDR3 recoveries.

Input: CDR3SeqData struct
Output: DataFrame

"""
function summarize(df::CDR3SeqData)
  statframe = @from i in df.counts begin
  @group i by i.patient into g
  @select {
      patient=first(g.patient),
      total_count=sum(g.count),
      min_count = minimum(g.count),
      max_count = maximum(g.count),
      avg_count = mean(g.count),
      median_count = median(g.count),
      num_unique_CDR3 = length(g.CDR3)
  }
  @collect DataFrame
  end

  return statframe

end

"""    filteredDataset(::Type{CDR3SeqData}, df::DataFrame, receptorname::String)

Subsets the dataset based on a receptor type of interest

Input: CDR3SeqData type, a DataFrame of counts, and a string that matches the types of the dataset
Output: DataFrame

Example: df_filtered = filteredDataset(CDR3SeqData, df, "TRA")

"""
function filteredDataset(::Type{CDR3SeqData}, df::DataFrame, receptorname::String)

  ## subset the type group only
  recepframe = @from i in df begin
    @where i.Receptor == receptorname
    @select i
    @collect DataFrame
  end;

  println("Number of recoveries: ", size(recepframe))
  println("Cohort size: ", length(unique(recepframe[:,:Filename])))

  ds = CDR3SeqData(recepframe)

  return ds
end

"""    filteredDataset(::Type{CDR3SeqData}, df::DataFrame, receptorname::Vector{String})

Subsets the dataset based on receptor types of interest, allowing users to subset based on combintations of types.

Input: CDR3SeqData type, a DataFrame of counts, and a string that matches the types of the dataset
Output: DataFrame

Example: df_filtered = filteredDataset(CDR3SeqData, df, ["TRA","TRB"])

"""
function filteredDataset(::Type{CDR3SeqData}, df::DataFrame, receptorname::Vector{String})
  ## subset the type group only
  recepframe = @from i in df begin
    @where i.Receptor in receptorname
    @select i
    @collect DataFrame
  end;

  println("Number of recoveries: ", size(recepframe))
  println("Cohort size: ", length(unique(recepframe[:,:Filename])))

  ds = CDR3SeqData(recepframe)

  return ds
end

"""    DiversityScores{T}

Stores a dataframe of diversity scores along with the range of qs that it was computed over.
"""
struct DiversityScores{T} <: Dataset
  qs::T
  data::DataFrame
end

DataFrame(ds::DiversityScores) = ds.data

"""    DiversityScores(countDF::DataFrame, rangeQ::AbstractVector, samplecol=:patient)

Computes point estimates of diversity from a countframe DataFrame, range of q values, and optional sample column identifier.

Point estimates of diversity of interest:
* low q diversity (q=0.01)
* high q diversity (q=100)
* delta qD diversity = low q diversity - high q diversity
* inflection point q
* magnitude of slope at the inflection point
* diversity qD evaluated at q=inflection point q
* Shannon diversity (log qD at q=1)
* Simpson diversity (1/qD at q=2)

Input: DataFrame of counts, vector of q values, optional sample column symbol (defaults to :patient)
Output: DiversityScores struct

"""
function DiversityScores(countDF::DataFrame, rangeQ::AbstractVector, samplecol=:patient)

  ## richness and evenness measures of diversity
  lowDiv = patientdiversity(countDF,rangeQ[1]);
  ids = lowDiv[:, samplecol];
  highDiv = patientdiversity(countDF,rangeQ[end]);

  ## inflection point related diversity metrics
  #    calcIPdiversity returns the q value of the inflection point & the slope at the inflection point
  dfIPslope = calcIPdiversity(countDF,rangeQ);
  # calculating the diversity qD at q=IP q value
  DatIPqDiv = Array{Float64}(undef,size(dfIPslope)[1])
  for row in 1:size(dfIPslope)[1]
    tmp = patientdiversity(countDF,dfIPslope[row,:InflPt]);
    DatIPqDiv[row,1] = tmp[row,:diversity]
  end

  ## commmonly used metrics of diversity
  ShanDiv = patientdiversity(countDF,1.0);
  SimpDiv = patientdiversity(countDF,2.0);

  ## compile all diversity scores into a single dataframe
  df = DataFrame(patient=ids,
     lowQ=lowDiv[:,:diversity],
     highQ=highDiv[:,:diversity],
     deltaqD=(lowDiv[:,:diversity].-highDiv[:,:diversity]),
     IPq=dfIPslope[:,:InflPt],
     IPslope=dfIPslope[:,:InflPtS],
     qDatIPq=DatIPqDiv,
     Shan=log.(ShanDiv[:,:diversity]),
     Simp= 1.0 ./(SimpDiv[:,:diversity]));
  return DiversityScores(rangeQ, df)
end

# function foo(c::SCRNASeqAnalysis, df::DataFrame)
#   #...SCRNASeq code here
# end

end

max_stars_repo_name: mcfefa/vdjDiversity

max_issues_repo_name: mcfefa/vdjDiversity

max_forks_repo_name: mcfefa/vdjDiversity


text:
# This is part of Calendars.jl. See the copyright note there.
# ======================== Islamic dates ====================

# The start of the Islamic calendar is AH 1, Muharram 1.
# That is equivalent to JD 622, July 16. We do not support
# dates prior to this date (no 'proleptic' Islamic calendar).
# Day number of the start of the Islamic calendar.
# AH-0001-01-01 = JD-0622-07-16 ~ DN#227015 = EpochIslamic
# AH-1111-11-11 = JD-1700-04-19 ~ DN#620667
# AH-6666-12-29 = JD-7089-11-28 ~ DN#2589222
const EpochIslamic = 227015
const ValidYearsIslamic = (1, 9666) # by convention

# True if year is an Islamic leap year.
function isLeapYearIslamic(year::DPart)
    return rem(11 * year + 14, 30) < 11
end

# Last month of Islamic year.
function LastMonthOfYearIslamic(year::DPart)
    isLeapYearIslamic(year) ? 13 : 12
end

# Last day in month during year on the Islamic calendar.
function LastDayOfMonthIslamic(year::DPart, month::DPart)
    b = (rem(month, 2) == 1
      || (month == 12 && isLeapYearIslamic(year)))
    return b ? 30 : 29
end

# Is the date a valid Islamic date?
function isValidDateIslamic(cd::CDate, warn=true)
    cal, year, month, day = cd
    println(cd)
    val = ( cal == AH
        && (ValidYearsIslamic[1] <= year && year <= ValidYearsIslamic[2])
        && (1 <= month && month <= LastMonthOfYearIslamic(year))
        && (1 <= day && day <= LastDayOfMonthIslamic(year, month)))
    ! val && warn && @warn(Warning(cd))
    return val
end

# Represents the day number a valid Islamic date?
function isValidDateIslamic(dn::DPart, warn=true)
    val = EpochIslamic <= dn
    ! val && warn && @warn(Warning(dn))
    return val
end

# Return the days this year so far.
# Does not depend on 'year' but kept in the signature.
function DayOfYearIslamic(year::DPart, month::DPart, day::DPart)
    # days so far this month plus ...
    day += ( 29 * (month - 1)   # days so far...
        + div(month, 2)         #   ...this year
    )
    return day
end

# Computes the day number from a valid Islamic date.
function DayNumberIslamic(year::DPart, month::DPart, day::DPart)

    day = DayOfYearIslamic(year, month, day)
    days = (EpochIslamic - 1      # days before start of calendar
         + 354 * (year - 1)       # non-leap days in prior years
         + div(3 + 11 * year, 30) # leap days in prior years
    )
    return day + days
end

# Computes the day number from a valid Islamic date.
DayNumberIslamic(d::CDate) = DayNumberIslamic(d[2], d[3], d[4])

# Computes the Islamic date from the day number.
function DateIslamic(dn::DPart)
    ! isValidDateIslamic(dn, false) && return InvalidDate

    # Search forward year by year from approximate year.
    year = div(dn - EpochIslamic, 355)
    while dn >= DayNumberIslamic(year + 1, 1, 1)
        year += 1
    end

    # Search forward month by month from Muharram.
    month = 1
    while dn > DayNumberIslamic(year, month,
                  LastDayOfMonthIslamic(year, month))
        month += 1
    end

    day = dn - DayNumberIslamic(year, month, 1) + 1

    return (AH, year, month, day)::CDate
end

# Return the day number of the first day in the given Islamic year.
function YearStartIslamic(year::DPart)
    return DayNumberIslamic(year, 1, 1)
end

# Return the day number of the last day in the given Islamic year.
function YearEndIslamic(year::DPart)
    return DayNumberIslamic(year + 1, 1, 1) - 1
end

# Return the number of months in the given Islamic year.
function MonthsInYearIslamic(year::DPart)
    return LastMonthOfYearIslamic(year)
end

# Return the number of days in the given Islamic year.
function DaysInYearIslamic(year::DPart)
    return YearStartIslamic(year + 1) - YearStartIslamic(year)
end

ext: jl


text:
config = MOI.Test.TestConfig()
# The test does not check the solution except the first variable so we just set zeros.
optimize!(mock) = MOIU.mock_optimize!(mock, [-0.184667; zeros(MOI.get(mock, MOI.NumberOfVariables()) - 1)])
for mock in mocks(optimize!)
    Tests.sos_quartic_comparison_test(mock, config)
end
optimize!(mock) = MOIU.mock_optimize!(mock, [-3.172412; zeros(MOI.get(mock, MOI.NumberOfVariables()) - 1)])
for mock in mocks(optimize!)
    Tests.sdsos_quartic_comparison_test(mock, config)
end
optimize!(mock) = MOIU.mock_optimize!(mock, [-11/3; zeros(MOI.get(mock, MOI.NumberOfVariables()) - 1)])
for mock in mocks(optimize!)
    Tests.dsos_quartic_comparison_test(mock, config)
end


text:
include("../NiLangCore/src/lens.jl")
using Lens

struct A
    a
    b::Float64
    c
    A() = new(0,0.0,0)
    A(x, y, z) = new(x, y, z)
end

a = A()
using BenchmarkTools
@benchmark @set $a.b = 3

@benchmark @with $a.b = 3
@benchmark chfield($a, Val(:b), 3)

@generated function chfield(x, ::Val{FIELD}, xval) where FIELD
    :(@with x.$FIELD = xval)
end


text:
# This approach heavily borrows from JSON3, with some changes
# for GraphQL specifics and other simplifications

# Each parse_{type} function
# - assumes that pos is the first character of what it is trying to parse (i.e. not an ignored charactaer)
# - returns pos at the position just after it has finished reading

"""
    parse(str::AbstractString)

Parses a GraphQL executable document string.
"""
function parse(str::AbstractString)
    buf = codeunits(str)
    len = length(buf)

    pos = 1
    line = 1
    column = 1
    @eof
    @skip_ignored

    definitions, pos, line, column = parse_definitions(buf, pos, line, column, len)
    return Document(definitions)
end

function parse_definitions(buf, pos, line, column, len)
    b = getbyte(buf, pos)
    definitions = Definition[]

    # Check for shorthand
    while pos <= len
        if b == UInt8('{')
            operation, pos, line, column = parse_operation(buf, pos, line, column, len)
            push!(definitions, operation)
        elseif !isnamestart(b)
            invalid("Expected name start character at start of definition", buf, pos)
        else
            # operation type
            definition_type, _, _, _ = parse_name(buf, pos, line, column, len)
            definition_type ∉ ("query", "mutation", "subscription", "fragment") && invalid("Definition type can't be \"$definition_type\"", buf, pos)
            if definition_type == "fragment"
                fragment, pos, line, column = parse_fragment_definition(buf, pos, line, column, len)
                push!(definitions, fragment)
            else
                operation, pos, line, column = parse_operation(buf, pos, line, column, len)
                push!(definitions, operation)
            end
        end
        only_ignored_left(buf, pos, len) && return definitions, pos, line, column
        @skip_ignored
    end
    return definitions, pos, line, column
end

function parse_operation(buf, pos, line, column, len)
    start_line = line
    start_column = column
    b = getbyte(buf, pos)
    if b == UInt8('{')
        operation_type = "query"
    elseif !isnamestart(b)
        invalid("Expected name start character", buf, pos)
    else
        # operation type
        operation_type, pos, line, column = parse_name(buf, pos, line, column, len)
        operation_type ∉ ("query", "mutation", "subscription") && throw(ArgumentError("operation type can't be $operation_type"))
    end

    @eof_skip_ignored

    if isnamestart(b)
        # name (optional)
        name, pos, line, column = parse_name(buf, pos, line, column, len)
    else
        name = nothing
    end

    @eof_skip_ignored

    if b == UInt('(')
        variables, pos, line, column = parse_variable_definitions(buf, pos, line, column, len)
    else
        variables = nothing
    end

    @eof_skip_ignored

    selection_set, pos, line, column = parse_selection_set(buf, pos, line, column, len)
    operation = Operation(operation_type, name, variables, nothing, selection_set, Loc(start_line, start_column))
    return operation, pos, line, column
end

function parse_fragment_definition(buf, pos, line, column, len)
    start_line = line
    start_column = column
    definition_type, pos, line, column = parse_name(buf, pos, line, column, len)
    definition_type != "fragment" && invalid("Fragment definition must start with 'fragment'", buf, pos)
    @eof_skip_ignored
    first_name_pos = pos
    name, pos, line, column = parse_name(buf, pos, line, column, len)
    @eof_skip_ignored
    on, pos, line, column = parse_name(buf, pos, line, column, len)
    on != "on" && invalid("Fragment definition must have form 'FragmentName on NamedTime'", buf, first_name_pos)
    @eof_skip_ignored
    named_type, pos, line, column = parse_name(buf, pos, line, column, len)
    @eof_skip_ignored
    if b == UInt('@')
        directives, pos, line, column = parse_directives(buf, pos, line, column, len)
    else
        directives = nothing
    end
    @eof_skip_ignored
    selection_set, pos, line, column = parse_selection_set(buf, pos, line, column, len)
    fragment_definition = FragmentDefinition(name, named_type, directives, selection_set, Loc(start_line, start_column))
    return fragment_definition, pos, line, column
end

function parse_variable_definitions(buf, pos, line, column, len)
    b = getbyte(buf, pos)
    if b != UInt('(') # redundant check?
        invalid("Variable definitions must start with '('", buf, pos)
    end

    pos += 1
    column += 1
    @eof_skip_ignored

    variables_defs = VariableDefinition[]
    while b != UInt(')')
        start_line = line
        start_column = column
        if b != UInt('$')
            invalid("Variable name must start with '\$'", buf, pos)
        end

        pos += 1
        column += 1
        @eof_skip_ignored

        name, pos, line, column = parse_name(buf, pos, line, column, len)

        @eof_skip_ignored

        if b != UInt(':')
            invalid("Variable name and type must be separated by ':'", buf, pos)
        end

        pos += 1  # Move past ':'
        column += 1
        @eof_skip_ignored

        type, pos, line, column = parse_type(buf, pos, line, column, len)

        @eof_skip_ignored

        # default value
        if b == UInt('=')
            pos += 1
            column += 1
            @eof_skip_ignored
            value, pos, line, column = parse_value(buf, pos, line, column, len)
        else
            value = nothing
        end

        # directives
        if b == UInt('@')
            pos += 1
            column += 1
            @eof_skip_ignored
            directive, pos, line, column = parse_directives(buf, pos, line, column, len)
        else
            directive = nothing
        end

        push!(variables_defs, VariableDefinition(name, type, value, directive, Loc(start_line, start_column)))
        @eof_skip_ignored
        b = getbyte(buf, pos)
    end
    pos += 1 # move past )
    column += 1
    return variables_defs, pos, line, column
end

function parse_directives(buf, pos, line, column, len)
    directives = Directive[]
    while pos <= len && getbyte(buf, pos) == UInt('@')
        directive, pos, line, column = parse_directive(buf, pos, line, column, len)
        push!(directives, directive)
        @eof_skip_ignored
    end
    return directives, pos, line, column
end

function parse_directive(buf, pos, line, column, len)
    start_line = line
    start_column = column
    b = getbyte(buf, pos)
    if b != UInt('@') # redundant check?
        invalid("Directives must start with '@'", buf, pos)
    end

    pos += 1
    column += 1
    @eof_skip_ignored

    name, pos, line, column = parse_name(buf, pos, line, column, len)

    if pos > len
        return Directive(name, nothing, Loc(start_line, start_column)), pos, line, column
    end

    @skip_ignored
    if b == UInt('(')
        arguments, pos, line, column = parse_arguments(buf, pos, line, column, len)
    end

    return Directive(name, arguments, Loc(start_line, start_column)), pos, line, column
end

function parse_type(buf, pos, line, column, len)
    type_str = ""
    b = getbyte(buf, pos)
    if b == UInt('[')
        type_str *= "["
        pos += 1
        column += 1
        @eof_skip_ignored
        name, pos, line, column = parse_type(buf, pos, line, column, len)
        type_str *= name
        @eof_skip_ignored
        b = getbyte(buf, pos)
        if b == UInt('!')
            type_str *= "!"
            pos += 1
            column += 1
            @eof_skip_ignored
        end
        b = getbyte(buf, pos)
        if b != UInt(']')
            invalid("Expected ']' in type", buf, pos)
        end
        type_str *= "]"
        pos += 1
        column += 1
    else
        name, pos, line, column = parse_name(buf, pos, line, column, len)
        type_str *= name
    end
    if pos > len
        return type_str, pos, line, column
    end
    b = getbyte(buf, pos)
    if b == UInt('!')
        type_str *= "!"
        pos += 1
        column += 1
    end
    return type_str, pos, line, column
end

function parse_name(buf, pos, line, column, len)
    b = getbyte(buf, pos)
    if !isnamestart(b)
        invalid("Expected name start character", buf, pos)
    end

    startpos = pos
    namelen = 0
    while isnamecontinue(b) && pos <= len
        pos += 1
        column += 1
        namelen += 1
        if pos <= len
            b = getbyte(buf, pos)
        end
    end
    name = unsafe_string(pointer(buf, startpos), namelen)
    return name, pos, line, column
end

function parse_selection_set(buf, pos, line, column, len)
    start_line = line
    start_column = column
    b = getbyte(buf, pos)
    @skip_ignored
    if b !== UInt8('{')
        invalid("Expected '{'", buf, pos)
    end

    pos += 1
    column += 1
    @eof_skip_ignored

    selections = Selection[]
    while b !== UInt8('}')
        selection, pos, line, column = parse_selection(buf, pos, line, column, len)
        push!(selections, selection)
        @eof
        @skip_ignored
    end
    pos += 1 # Move past '}'
    column += 1

    return SelectionSet(selections, Loc(start_line, start_column)), pos, line, column
end

function parse_selection(buf, pos, line, column, len)
    b = getbyte(buf, pos)
    if b == UInt('.') && getbyte(buf, pos+1) == UInt('.') && getbyte(buf, pos+2) == UInt('.')
        return parse_fragment(buf, pos, line, column, len)
    else
        # Field
        return parse_field(buf, pos, line, column, len)
    end
end

function parse_fragment(buf, pos, line, column, len)
    b = getbyte(buf, pos)
    if !(b == UInt('.') && getbyte(buf, pos+1) == UInt('.') && getbyte(buf, pos+2) == UInt('.'))
        # Are these and similar conditions necessary if we're already checking further up?
        invalid("Inline fragments and fragment spread must start '...", buf, pos)
    end
    pos += 3
    column += 3
    start_line = line
    start_column = column
    @eof_skip_ignored
    b = getbyte(buf, pos)
    if isnamestart(b)
        name, pos, line, column = parse_name(buf, pos, line, column, len)
        if name !== "on"
            @eof_skip_ignored
            b = getbyte(buf, pos)
            if b == UInt('@')
                directives, pos, line, column = parse_directives(buf, pos, line, column, len)
            else
                directives = nothing
            end
            @eof_skip_ignored
            b = getbyte(buf, pos)
            return FragmentSpread(name, directives, Loc(start_line, start_column)), pos, line, column
        end

        # Inline fragment starting with type condition
        @eof_skip_ignored
        named_type, pos, line, column = parse_name(buf, pos, line, column, len)
        @eof_skip_ignored
    else
        named_type = nothing
    end

    # Only get here with inline fragment
    if b == UInt('@')
        directives, pos, line, column = parse_directives(buf, pos, line, column, len)
    else
        directives = nothing
    end

    @eof_skip_ignored

    selection_set, pos, line, column = parse_selection_set(buf, pos, line, column, len)
    return InlineFragment(named_type, directives, selection_set, Loc(start_line, start_column)), pos, line, column
end

function parse_number(buf, pos, line, column, len)
    b = getbyte(buf, pos)

    number_start = pos
    number_length = 0

    if b == UInt('-')
        pos += 1
        column += 1
        number_length += 1
        @eof
        b = getbyte(buf, pos)
    end

    if b == UInt('0')
        if pos + 1 < len && isdigit(Char(getbyte(buf, pos+1)))
            error("IntValue or FloatValue cannot begin with leading zero")
        end
    end

    while isdigit(Char(b))
        number_length += 1
        pos += 1
        column += 1
        if pos > len
            # Need this check when single digit is at end of buffer
            @inbounds @views value = Parsers.parse(Int64, buf[number_start:number_start+number_length-1])
            return value, pos, line, column
        end
        b = getbyte(buf, pos)
    end

    if !(b == UInt('.') || isexponent(b))
        if islowercaseletter(b) || iscapitalletter(b)
            error("Integer value cannot be followed by a name char") # Refine
        end
        @inbounds @views value = Parsers.parse(Int64, buf[number_start:number_start+number_length-1])
    else
        while isdigit(Char(b)) || b == UInt('.') || isexponent(b) || b == UInt('+') || b == UInt('-')
            # Could do better checking here
            number_length += 1
            pos += 1
            column += 1
            if pos > len
                # shouldn't get here but useful for testing function
                break
            end
            b = getbyte(buf, pos)
        end
        @inbounds @views value = Parsers.parse(Float64, buf[number_start:number_start+number_length-1])
    end
    return value, pos, line, column
end

function parse_string(buf, pos, line, column, len)
    b = getbyte(buf, pos)
    escaped = false
    if b == UInt('"') && getbyte(buf, pos+1) == UInt('"') && getbyte(buf, pos+2) == UInt('"')
        isblockstr = true
        pos += 3
        column += 3
        @eof
        b = getbyte(buf, pos)
        stringlen = 0
        stringstart = pos

        while !(b == UInt('"') && getbyte(buf, pos+1) == UInt('"') && getbyte(buf, pos+2) == UInt('"'))
            if b == UInt('\\') && all((i) -> getbyte(buf, pos+i) == UInt('"'), (1,2,3))
                # Escaped block string
                escaped = true
                pos += 4
                stringlen += 4
                column += 4
            else
                if islineterminator(b)
                    line += 1
                    column = 1
                else
                    column += 1
                end
                pos += 1
                stringlen += 1
            end
            @eof
            b = getbyte(buf, pos)
        end
        # Get past remaining """
        pos += 3
        column += 3
    else
        isblockstr = false
        pos += 1
        column += 1
        @eof
        b = getbyte(buf, pos)
        stringlen = 0
        stringstart = pos
        while b != UInt('"')
            if b == UInt('\\')
                escaped = true
                if getbyte(buf, pos+1) == UInt('u')
                    # unicode
                    if !all((i) -> isunicodechar(getbyte(buf, pos+i)), (2,3,4,5))
                        invalid("Unicode must have four digits", buf, pos)
                    end
                    pos += 5
                    column += 5
                    stringlen += 5
                else
                    # other escaped characted
                    pos += 2
                    column += 2
                    stringlen += 2
                end
            else
                pos += 1
                column += 1
                stringlen += 1
            end
            @eof
            b = getbyte(buf, pos)
        end

        # Get past remaining "
        pos += 1
        column += 1
    end
    if escaped
        if isblockstr
            return format_block_string(unescape_blockstr(PointerString(pointer(buf, stringstart), stringlen))), pos, line, column
        else
            return unescape(PointerString(pointer(buf, stringstart), stringlen)), pos, line, column
        end
    else
        if isblockstr
            return format_block_string(unsafe_string(pointer(buf, stringstart), stringlen)), pos, line, column
        else
            return unsafe_string(pointer(buf, stringstart), stringlen), pos, line, column
        end
    end
end

function parse_arguments(buf, pos, line, column, len)
    b = getbyte(buf, pos)
    if b != UInt('(')
        invalid("Arguments must start with (", buf, pos)
    end

    pos += 1
    column += 1
    @eof_skip_ignored

    arguments = Argument[]

    while b!= UInt(')')
        argument, pos, line, column = parse_argument(buf, pos, line, column, len)
        push!(arguments, argument)
        @eof_skip_ignored
    end

    isempty(arguments) && invalid("Expected at least one argument", buf, pos)
    # Move past ')'
    pos +=1
    column += 1

    return arguments, pos, line, column
end
function parse_argument(buf, pos, line, column, len)
    start_line = line
    start_column = column
    # Get name (either name)
    name, pos, line, column = parse_name(buf, pos, line, column, len)
    b = getbyte(buf, pos)
    @skip_ignored

    # Check :
    if b != UInt(':')
        invalid("Expected : in argument", buf, pos)
    end
    pos += 1
    column += 1
    @eof_skip_ignored

    # Get value
    value, pos, line, column = parse_value(buf, pos, line, column, len)

    return Argument(name, value, Loc(start_line, start_column)), pos, line, column
end

function parse_value(buf, pos, line, column, len)
    b = getbyte(buf, pos)
    if isdigit(Char(b)) || b == UInt('-')
        value, pos, line, column = parse_number(buf, pos, line, column, len)
    elseif b == UInt('t') && getbyte(buf, pos+1) == UInt('r') && getbyte(buf, pos+2) == UInt('u') && getbyte(buf, pos+3) == UInt('e')
        # will fail if shorter than this
        value = true
        pos += 4
        column += 4
    elseif b == UInt('f') && getbyte(buf, pos+1) == UInt('a') && getbyte(buf, pos+2) == UInt('l') && getbyte(buf, pos+3) == UInt('s') && getbyte(buf, pos+4) == UInt('e')
        value = false
        pos += 5
        column += 5
    elseif b == UInt('n') && getbyte(buf, pos+1) == UInt('u') && getbyte(buf, pos+2) == UInt('l') && getbyte(buf, pos+3) == UInt('l')
        value = nothing
        pos += 4
        column += 4
    elseif b == UInt('"')
        value, pos, line, column = parse_string(buf, pos, line, column, len)
        @eof # needed as not checked in parse_string
    elseif b == UInt('[')
        value, pos, line, column = parse_list(buf, pos, line, column, len)
    elseif b == UInt('{')
        # Input object
        value, pos, line, column = parse_input_object(buf, pos, line, column, len)
    elseif isnamestart(b)
        # enum
        start_line = line
        start_column = column
        value_str, pos, line, column = parse_name(buf, pos, line, column, len)
        value = Enum(value_str, Loc(start_line, start_column))
    elseif b == UInt('$')
        # variable
        variable_location = Loc(line, column)
        pos += 1
        column += 1
        @eof_skip_ignored
        name, pos, line, column = parse_name(buf, pos, line, column, len)
        value = Variable(name, variable_location)
    else
        invalid("Could not parse value", buf, pos)
    end

    return value, pos, line, column
end

function parse_input_object(buf, pos, line, column, len)
    start_line = line
    start_column = column
    b = getbyte(buf, pos)
    if b != UInt('{')
        invalid("Input object must start with {", buf, pos)
    end

    pos += 1
    column += 1
    @eof_skip_ignored

    object_fields = ObjectField[]

    while b != UInt('}')
        field_start_line = line
        field_start_column = column
        if !isnamestart(b)
            invalid("Expected name start character or '}'", buf, pos)
        end

        name, pos, line, column = parse_name(buf, pos, line, column, len)

        @eof_skip_ignored

        if b != UInt(':')
            invalid("ObjectField name must be followed by ':'", buf, pos)
        end

        pos += 1
        column += 1
        @eof_skip_ignored

        value, pos, line, column = parse_value(buf, pos, line, column, len)
        push!(object_fields, ObjectField(name, value, Loc(field_start_line, field_start_column)))
        @eof_skip_ignored
    end

    # Move past '}'
    pos += 1
    column += 1

    return InputObject(object_fields, Loc(start_line, start_column)), pos, line, column
end

function parse_list(buf, pos, line, column, len)
    b = getbyte(buf, pos)
    if b != UInt('[')
        invalid("List must start with [", buf, pos)
    end

    pos += 1
    column += 1
    @eof_skip_ignored
    b = getbyte(buf, pos)

    list = Any[] # TODO: we'd use introspection here to determine the type of the list.

    first_element = true
    while b != UInt(']')
        value, pos, line, column = parse_value(buf, pos, line, column, len)
        push!(list, value)
        !first_element && @assert typeof(value) == typeof(first(list))
        @eof_skip_ignored
    end

    # Move past ']'
    pos += 1
    column += 1

    return list, pos, line, column
end

function parse_field(buf, pos, line, column, len)
    start_line = line
    start_column = column
    # Get first name (either name or alias)
    name, pos, line, column = parse_name(buf, pos, line, column, len)
    alias::Union{String, Nothing} = nothing
    if pos > len
        # Shouldn't really happen but useful for tests
        field = Field(alias, name, nothing, nothing, nothing, Loc(start_line, start_column))
        return field, pos, line, column, len
    end
    b = getbyte(buf, pos)
    @skip_ignored

    # Check for alias
    if b == UInt(':')
        # Alias
        pos += 1
        column += 1
        @eof_skip_ignored
        alias = name
        name, pos, line, column = parse_name(buf, pos, line, column, len)
    end
    if pos > len
        # Shouldn't really happen but useful for tests
        field = Field(alias, name, nothing, nothing, nothing, Loc(start_line, start_column))
        return field, pos
    end
    b = getbyte(buf, pos)
    @skip_ignored

    # Check for arguments
    if b == UInt8('(')
        arguments, pos, line, column = parse_arguments(buf, pos, line, column, len)
    else
        arguments = nothing
    end
    b = getbyte(buf, pos)
    @skip_ignored # TODO: what if arguments and nothing else after?

    # Check for directives
    if b == UInt('@')
        directives, pos, line, column = parse_directives(buf, pos, line, column, len)
    else
        directives = nothing
    end
    b = getbyte(buf, pos)
    @skip_ignored # TODO: what if arguments and nothing else after?

    # Check for selection set
    if b == UInt8('{')
        # SelectionSet
        selection_set, pos, line, column = parse_selection_set(buf, pos, line, column, len)
    else
        selection_set = nothing
    end

    # Build operation from tape
    field = Field(alias, name, arguments, directives, selection_set, Loc(start_line, start_column))
    return field, pos, line, column
end


text:
import LambdaCalculus
using LambdaCalculus.DeBruijn
using PEG

export encode,
    decode,
    compress,
    decompress

encode(var::Var) = "\x01" ^ var.index * "\x00"
encode(abs::Abs) = "\x00\x00" * encode(abs.body)
encode(app::App) = "\x00\x01" * encode(app.car) * encode(app.cdr)

@rule term = var, abs, app
@rule var = "\x01"[+] & "\x00" > ((i, _) -> Var(length(i)))
@rule abs = "\x00\x00" & term > ((_, y) -> Abs(y))
@rule app = "\x00\x01" & term & term > ((_, car, cdr) -> App(car, cdr))

"""
    decode(input::AbstractString)

Parse a De Bruijn indexed term from its binary representation.
"""
function decode(input::AbstractString; onsuccess = identity, onerror = throw)
    try
        onsuccess(parse_whole(term, input))
    catch parser_error
        onerror(parser_error)
    end
end

compress(term::Term) = BitVector([c == '\x01' for c in encode(term)])

function decompress(compressed::BitVector; onsuccess = identity, onerror = throw)
    decode(String(map(Char, compressed)); onsuccess = onsuccess, onerror = onerror)
end


text:
import Base: hcat, vcat, \
export hcat, vcat, \, diagOp, A_mul_B

include("Composition.jl")
include("NFFTOp.jl")
include("ExplicitOp.jl")
include("SensitivityOp.jl")
include("SamplingOp.jl")
include("MapSliceOp.jl")
include("FieldmapNFFTOp.jl")
include("EncodingOp.jl")
include("SparseOp.jl")

"""
    hcat(A::AbstractLinearOperator, n::Int)

horizontally concatenates alinear operator `A` n times
"""
function hcat(A::AbstractLinearOperator, n::Int)
  op = A
  for i = 2:n
    op = LinearOperators.hcat(op, A)
  end
  return op
end

"""
    hcat(A::AbstractLinearOperator, n::Int)

vertically concatenates alinear operator `A` n times
"""
function vcat(A::AbstractLinearOperator, n::Int)
  op = A
  for i = 2:n
    op = LinearOperators.vcat(op, A)
  end
  return op
end

function diagOpProd(x::Vector{T}, nrow::Int, xIdx, yIdx, ops :: AbstractLinearOperator...) where T
  y = Vector{T}(undef, nrow)
  @sync for i=1:length(ops)
    Threads.@spawn begin
      y[yIdx[i]:yIdx[i+1]-1] = ops[i]*x[xIdx[i]:xIdx[i+1]-1]
    end
  end
  return y
end

function diagOpTProd(x::Vector{T}, ncol::Int, xIdx, yIdx, ops :: AbstractLinearOperator...) where T
  y = Vector{T}(undef, ncol)
  @sync for i=1:length(ops)
    Threads.@spawn begin
      y[yIdx[i]:yIdx[i+1]-1] = transpose(ops[i])*x[xIdx[i]:xIdx[i+1]-1]
    end
  end
  return y
end

function diagOpCTProd(x::Vector{T}, ncol::Int, xIdx, yIdx, ops :: AbstractLinearOperator...) where T
  y = Vector{T}(undef, ncol)
  @sync for i=1:length(ops)
    Threads.@spawn begin
      y[yIdx[i]:yIdx[i+1]-1] = adjoint(ops[i])*x[xIdx[i]:xIdx[i+1]-1]
    end
  end
  return y
end

function Base.copy(S::LinearOperator{T}) where T
  deepcopy(S)
end

mutable struct DiagOp{T} <: AbstractLinearOperator{T}
  nrow :: Int
  ncol :: Int
  symmetric :: Bool
  hermitian :: Bool
  prod! :: Function
  tprod! :: Function
  ctprod! :: Function
  nprod :: Int
  ntprod :: Int
  nctprod :: Int
  args5 :: Bool
  use_prod5! :: Bool
  allocated5 :: Bool
  Mv5 :: Vector{T}
  Mtu5 :: Vector{T}
  ops
  equalOps :: Bool
  xIdx :: Vector{Int}
  yIdx :: Vector{Int}
end

"""
    diagOp(ops :: AbstractLinearOperator...)

create a bloc-diagonal operator out of the `LinearOperator`s contained in ops
"""
function diagOp(ops :: AbstractLinearOperator...)
  nrow = 0
  ncol = 0
  S = eltype(ops[1])
  for i = 1:length(ops)
    nrow += ops[i].nrow
    ncol += ops[i].ncol
    S = promote_type(S, eltype(ops[i]))
  end

  xIdx = cumsum(vcat(1,[ops[i].ncol for i=1:length(ops)]))
  yIdx = cumsum(vcat(1,[ops[i].nrow for i=1:length(ops)]))

  Op = DiagOp{S}( nrow, ncol, false, false,
                     (res,x) -> (res .= diagOpProd(x,nrow,xIdx,yIdx,ops...)),
                     (res,y) -> (res .= diagOpTProd(y,ncol,yIdx,xIdx,ops...)),
                     (res,y) -> (res .= diagOpCTProd(y,ncol,yIdx,xIdx,ops...)),
                     0, 0, 0, false, false, false, S[], S[],
                     [ops...], false, xIdx, yIdx)

  return Op
end

function diagOp(op::AbstractLinearOperator, N=1)
  nrow = N*op.nrow
  ncol = N*op.ncol
  S = eltype(op)
  ops = [copy(op) for n=1:N]

  xIdx = cumsum(vcat(1,[ops[i].ncol for i=1:length(ops)]))
  yIdx = cumsum(vcat(1,[ops[i].nrow for i=1:length(ops)]))

  Op = DiagOp{S}( nrow, ncol, false, false,
                    (res,x) -> (res .= diagOpProd(x,nrow,xIdx,yIdx,ops...)),
                    (res,y) -> (res .= diagOpTProd(y,ncol,yIdx,xIdx,ops...)),
                    (res,y) -> (res .= diagOpCTProd(y,ncol,yIdx,xIdx,ops...)),
                     0, 0, 0, false, false, false, S[], S[],
                     ops, true, xIdx, yIdx )

  return Op
end

### Normal Matrix Code ###

struct DiagNormalOp{U,V,T}
  ops::U
  normalOps::V
  nrow::Int64
  ncol::Int64
  idx::Vector{Int64}
  y::Vector{T}
end

function SparsityOperators.normalOperator(S::DiagOp, W=I)
  weights = W*ones(S.nrow)

  T = eltype(S)

  if S.equalOps
    # this opimization is only allowed if all ops are the same
    opInner = normalOperator(S.ops[1], WeightingOp(weights[S.yIdx[1]:S.yIdx[2]-1].^2))
    op = DiagNormalOp(S.ops, [copy(opInner) for i=1:length(S.ops)], S.ncol, S.ncol, S.xIdx, zeros(T, S.ncol) )
  else
    op = DiagNormalOp(S.ops, [normalOperator(S.ops[i], WeightingOp(weights[S.yIdx[i]:S.yIdx[i+1]-1].^2))
                     for i in 1:length(S.ops)], S.ncol, S.ncol, S.xIdx, zeros(T, S.ncol) )
  end

  return op
end

function Base.:*(S::DiagNormalOp, x::AbstractVector{T}) where T
  _produ_diagnormalop(S.normalOps, S.idx, x, S.y)
  return S.y
end

function _produ_diagnormalop(ops, idx, x::AbstractVector{T}, y::AbstractVector{T}) where T
  @sync for i=1:length(ops)
    Threads.@spawn begin
       y[idx[i]:idx[i+1]-1] = ops[i]*x[idx[i]:idx[i+1]-1]
    end
  end
  return
end

#
# use hermitian conjugate as an estimate for the inverse (fallback)
#
\(A::AbstractLinearOperator, x::Vector) = adjoint(A)*x


text:
#!/usr/bin/julia
# FizzBuzz using basic division rules
# Author: @joaogui1

function div3(n)
    s = string(n)
    while length(s) > 1
        s = string(sum([parse(Int, c) for c in s]))
    end
    return s in ["0", "3", "6", "9"]
end

function div5(n)
    s = string(n)
    return s[end] in "05"
end

function fizzbuzz(n)
    for i = 1:n
        print(i, ": ")
        by3, by5 = div3(i), div5(i)

        if by3 && by5
            println("FizzBuzz")
        elseif by3
            println("Fizz")
        elseif by5
            println("Buzz")
        else
            println(i)
        end
    end
end

fizzbuzz(100)


text:
export  SQLAnd, SQLOr, SQLNot

processentry( x ) = "$x"
processentry( x::Vector ) = string( '[', join( processentry.(x), ", " ), ']' )
processentry( x::Missing ) = "NULL"
processentry( x::Dict ) = json(x)

sqlcops = Dict(
    :lt => (2, ["<"], :geq),
    :gt => (2, [">"], :leq),
    :leq => (2, ["<="], :gt),
    :geq => (2, [">="], :lt),
    :eq => (2, ["="], :neq),
    :neq => (2, ["<>"], :eq),
    :in => (2, ["IN"], :nin),
    :nin => (2, ["NOT IN"], :in),
    :like => (2, ["LIKE"], :nlike),
    :nlike => (2, ["NOT LIKE"], :like),
    :between => (3, ["BETWEEN", "AND"], :nbetween),
    :nbetween => (3, ["NOT BETWEEN", "AND"], :between),
    :betweensym => (3, ["BETWEEN SYMMETRIC", "AND"], :nbetweensym),
    :nbetweensym => (3, ["NOT BETWEEN SYMMETRIC", "AND"], :betweensym),
    :distinct => (2, ["IS DISTINCT FROM"], :ndistinct),
    :ndistinct => (2, ["IS NOT DISTINCT FROM"], :distinct),
    :null => (1, ["IS NULL"], :nnull),
    :nnull => (1, ["IS NOT NULL"], :null),
    :true => (1, ["IS TRUE"], :ntrue),
    :ntrue => (1, ["IS NOT TRUE"], :true),
    :false => (1, ["IS FALSE"], :nfalse),
    :nfalse => (1, ["IS NOT FALSE"], :false),
    :unknown => (1, ["IS UNKNOWN"], :nunknown),
    :nunknown => (1, ["IS NOT UNKNOWN"], :unknown)
)

abstract type SQLCondition end

struct NoCondition <: SQLCondition end
"""
`nocond` is the empty SQL condition. This will not render if part of a clause.
"""
const nocond = NoCondition()
Base.isempty(cond::SQLCondition) = cond isa NoCondition

mutable struct CompoundCondition <: SQLCondition
    lop::Symbol
    conds::Vector{SQLCondition}

    function CompoundCondition( lop, conds::Vector{T} ) where T <: SQLCondition
        lop ∈ [:and, :or, :not] || @error( "Unknown logical operator, only AND, OR, and NOT allowed." )
        length(conds) == 0 && @error( "At least one SQL condition must be given." )

        if lop === :not
            length(conds) == 1 || @error( "Exactly one SQL must be given with NOT operator." )
            return new( :not, deepcopy(conds) )
        end

        length(conds) == 1 && return conds[1]
        new( lop, deepcopy(conds) )
    end

    CompoundCondition( lop::Symbol, conds::SQLCondition... ) = CompoundCondition( lop, conds |> collect )
end

function Base.show( io::IO, cc::CompoundCondition )
    if cc.lop === :not
        print( io, "NOT ($(cc.conds[1]))" )
        return
    end

    cstrs = map( cc.conds ) do cond
        cond isa CompoundCondition ? "($cond)" : "$cond"
    end
    print( io, join( cstrs, " $(string(cc.lop) |> uppercase) " ))
end

"""
```
SQLAnd( conds::Vector{T} ) where T <: SQLCondition
SQLAnd( conds::SQLCondition... )
```
This function creates an SQL condition that concatenates the SQL conditions in `conds` with the `AND` operator.
"""
SQLAnd( conds::Vector{T} ) where T <: SQLCondition = CompoundCondition( :and, conds )
SQLAnd( conds::SQLCondition... ) = CompoundCondition( :and, conds... )
"""
```
SQLOr( conds::Vector{T} ) where T <: SQLCondition
SQLOr( conds::SQLCondition... )
```
This function creates an SQL condition that concatenates the SQL conditions in `conds` with the `OR` operator.
"""
SQLOr( conds::Vector{T} ) where T <: SQLCondition = CompoundCondition( :or, conds )
SQLOr( conds::SQLCondition... ) = CompoundCondition( :or, conds... )
"""
```
SQLNot( cond::CompoundCondition )
SQLNot( cond::BaseCondition )
```
This function creates an SQL condition that is the logical `NOT` of the SQL condition in `cond`. If a basic SQL condition is passed as argument, it creates a new basic condition with the inverse of the operator.
"""
SQLNot( cond::CompoundCondition ) = CompoundCondition( :not, [cond] )

mutable struct BaseCondition <: SQLCondition
    cop::Symbol
    args::Vector{String}

    function BaseCondition( cop::Symbol, args::Vector )
        haskey( sqlcops, cop ) || @error "Unknown comparison operator."
        nargs = length(args)
        rargs = sqlcops[cop][1]
        if cop ∈ [:in, :nin]
            nargs >= rargs || @error "Incorrect number of arguments, at least $rargs arguments required."
        else
            nargs == rargs || @error "Incorrect number of arguments, $rargs arguments required."
        end
        new( cop, processentry.(args) )
    end

    BaseCondition( cop::Symbol, args... ) = BaseCondition( cop, args |> collect )
end

function Base.show( io::IO, bc::BaseCondition )
    nargs = length(bc.args)
    cparts = fill( "", nargs == 1 ? 2 : 2*nargs - 1 )
    cparts[1:2:end] = bc.args
    cparts[2:2:end] = sqlcops[bc.cop][2]
    bc.cop ∈ [:in, :nin] && (cparts[3] = amendarg( cparts[3] ))
    print( io, join( cparts, " " ) )
end

amendarg( arg::String ) =  string( '(', arg[(endswith( arg, ";" ) ? 1 : 2):(end-1)], ')' )

SQLNot( cond::BaseCondition ) = BaseCondition( sqlcops[cond.cop][3], cond.args... )

for cop in filter( cop -> sqlcops[cop][1] == 1, keys(sqlcops) |> collect )
    etype = cop ∈ [:null, :nnull] ? Union{Real, SVar, SQLString} : Union{Bool, SVar, SQLCondition, SQLString}
    Core.eval( @__MODULE__, "export SQL$cop" |> Meta.parse )
    Core.eval( @__MODULE__, """\"\"\"
    ```
    SQL$cop( arg::$etype )
    ```
    This function creates a SQL condition of type `<arg> $(sqlcops[cop][2][1])`.
    \"\"\"
    SQL$cop( arg::$etype ) = BaseCondition( :$cop, arg )""" |> Meta.parse )
end

for cop in filter( cop -> sqlcops[cop][1] == 2, keys(sqlcops) |> collect )
    etype = Union{Real, SVar, SQLString}
    Core.eval( @__MODULE__, "export SQL$cop" |> Meta.parse )

    if cop ∈ [:in, :nin]
        Core.eval( @__MODULE__, """\"\"\"
        ```
        SQL$cop(
            larg::$etype,
            rarg::Vector{T}
        ) where T <: $etype
        SQL$cop(
            larg::$etype,
            rarg::SQLStatement )
        ```
        This function creates a SQL condition of type `<larg> $(sqlcops[cop][2][1]) <rarg>`.
        \"\"\"
        SQL$cop( larg::$etype, rarg::Vector{T} ) where T <: $etype = BaseCondition( :$cop, larg, rarg )""" |> Meta.parse )
        Core.eval( @__MODULE__, """SQL$cop( larg::$etype, rarg::SQLStatement ) = BaseCondition( :$cop, larg, rarg )""" |> Meta.parse )
    elseif cop ∈ [:like, :nlike]
        Core.eval( @__MODULE__, """\"\"\"
        ```
        SQL$cop(
            larg::$etype,
            rarg::Union{AbstractString, SQLString} )
        ```
        This function creates a SQL condition of type `<larg> $(sqlcops[cop][2][1]) <rarg>`.
        \"\"\"
        SQL$cop( larg::$etype, rarg::Union{AbstractString, SQLString} ) = BaseCondition( :$cop, larg, rarg |> SQLString )""" |> Meta.parse )
    else
        Core.eval( @__MODULE__, """\"\"\"
        ```
        SQL$cop(
            larg::$etype,
            rarg::$etype )
        ```
        This function creates a SQL condition of type `<larg> $(sqlcops[cop][2][1]) <rarg>`.
        \"\"\"
        SQL$cop( larg::$etype, rarg::$etype ) = BaseCondition( :$cop, larg, rarg )""" |> Meta.parse )
    end
end

for cop in filter( cop -> sqlcops[cop][1] == 3, keys(sqlcops) |> collect )
    etype = Union{Real, SVar, SQLString}
    Core.eval( @__MODULE__, "export SQL$cop" |> Meta.parse )
    Core.eval( @__MODULE__, """\"\"\"
    ```
    SQL$cop(
        carg::$etype,
        larg::$etype,
        rarg::$etype )
    ```
    This function creates a SQL condition of type `<carg> $(sqlcops[cop][2][1]) <larg> $(sqlcops[cop][2][2]) <rarg>`.
    \"\"\"
    SQL$cop( carg::$etype, larg::$etype, rarg::$etype ) = BaseCondition( :$cop, carg, larg, rarg )""" |> Meta.parse )
end


text:
function FeatureCollection(url)
    ch = Channel{STAC.Item}() do c
        while true
            @debug "get $url"
            r = HTTP.get(url)
            data = JSON3.read(String(r.body))
            for d in data[:features]
                put!(c,STAC.Item("",d,STAC._assets(d)))
            end

            # check if there is a next page
            next = filter(d -> get(d,"rel",nothing) == "next",data[:links])
            # no more next page
            if length(next) == 0
                break
            else
                url = next[1][:href]
            end
        end
    end
end

"""
    search(cat::Catalog, collections, lon_range, lat_range, datetime; limit = 200)

Search items in a `STAC.Catalog` `cat` belong to the `collections` within the
longitude range `lon_range` (start and end longitude),
latitude range `lat_range`, and time range (start and end `DateTime`).

The function `search` returns an julia `Channel` with the search results over
which one can only iterate once. If the results should be saved use
`search_results = collect(search(cat,...))`.

Example:

```julia
using STAC, Dates
collections = "landsat-8-c2-l2"
time_range = (DateTime(2018,01,01), DateTime(2018,01,02))
lon_range = (2.51357303225, 6.15665815596)
lat_range = (49.5294835476, 51.4750237087)

catalog = STAC.Catalog("https://planetarycomputer.microsoft.com/api/stac/v1")

search_results = collect(search(catalog, collections, lon_range, lat_range, time_range))
```

"""
function search(cat::Catalog, collections, lon_range, lat_range, time_range; limit = 200)
    format_datetime(x::Union{Dates.Date,Dates.DateTime}) = Dates.format(x,Dates.dateformat"yyyy-mm-ddTHH:MM:SS.sssZ")
    format_datetime(x::Union{AbstractVector,Tuple}) = join(format_datetime.(x),"/")

    west, east = lon_range
    south, north = lat_range
    bbox = join(string.([west, south, east, north]),",")
    query = Dict(
        "collections" => collections,
        "bbox" => bbox,
        "limit" => limit,
        "datetime" => format_datetime(time_range),
    )

    @debug "query: $query"
    url = string(URI(URI(cat.url * "/search"), query = query))

    return FeatureCollection(url)
end

export search

max_stars_repo_name: Alexander-Barth/STACatalog

max_issues_repo_name: Alexander-Barth/STACatalog

max_forks_repo_name: Alexander-Barth/STACatalog


text:
break_on(:error)

"""
    teralasso(S_kGram, regtype, regparam, regcoef, tol, niter, ninner)

Implementation of Tensor Graphical Lasso (TeraLasso) precision matrix estimator.
Learns a K-order Kronecker sum covariance model.

# Arguments
- `S_kGram::AbstractVector`: length K array of factorwise
    Gram matrices (S_k in the paper).
    NOTE: use the most general accepted type (UnionAll):
    ```julia-repl
    julia> AbstractVector{<:AbstractMatrix{<:Real}} <: AbstractVector === AbstractArray{T, 1} where T
    true

    julia> AbstractVector{<:AbstractMatrix{<:Real}} <: AbstractVector{Any}
    false
    ```
- `regtype::Symbol=:L1`: type of regularization, :L1, :SCAD, or :MCP.
- `regparam::Real=20`: parameter for SCAD/MCP (a in the paper).
- `regcoef::AbstractVector`: length K vector of regularization
    coefficients for each factor (λ).

# Keywords
- `zeta::Real=0.1`: initial step size.
- `c::Real=0.5`: backtracking constant.
    NOTE: 0.5 too large for active region data; use 0.05.
- `tol::Real=1e-6`: tolerance for convergence criterion (ε).
- `niter::Int=50`: maximum number of iterations allowed.
- `ninner::Int=1`: maximum number of backtracking line search iteration.

# Returns:
- `Psi::AbstractVector{AbstractMatrix{Real}}`: length K array of factorwise
    precision matrices
- `out::AbstractVector{AbstractVector{Any}}`: iteration history.
"""

function teralasso(
    S_kGram::AbstractVector,
    regtype::Symbol=:L1,
    regparam::Real=20,
    regcoef::Union{Nothing, AbstractVector}=nothing,
    ;
    zeta::Real=0.1,
    c::Real=0.1,
    tol::Real=1e-6,
    niter::Int=50,
    ninner::Int=11,
)
    # Set eigenvalue upper bound for nonconvex contraints.
    if regtype == :L1
        mu = 1e-6
    elseif regtype == :SCAD
        mu = 1 / (regparam - 1)
    elseif regtype == :MCP
        mu = 1 / regparam
    end
    kappa = sqrt(2 / mu) * 2

    # Preparation
    K = length(S_kGram)
    ps = [size(S_kGram[k], 1) for k in 1:K]
    p = prod(ps)
    regcoef = isnothing(regcoef) ? sqrt.(log.(ps) ./ (p ./ ps)) : regcoef
    #regcoef .+= 1e-9  # inv(cov) may be ill conditioned if p < n and lambda = 0
    S_tilde = [S - (K - 1) / K * mean(diag(S)) * I for S in S_kGram]

    # Initlalization
    logdet = 0
    logdetN = 0 # updated outside inner loop
    Psi = 1 / K * [Matrix{Float64}(I, ps[k], ps[k]) for k = 1:K]
    PsiN = 1 / K * [Matrix{Float64}(I, ps[k], ps[k]) for k = 1:K]
    M = 1 / K * [Matrix{Float64}(I, ps[k], ps[k]) for k = 1:K]
    out = Array{Any}(undef, niter)

    for iter = 1:niter
        # Backtracking line search
        for cnt = 1:ninner
            for k = 1:K
                PsiN[k] = shrink_regularizer(
                    Psi[k] - zeta * (S_tilde[k] - M[k]),
                    regcoef[k],
                    zeta,
                    regtype,
                    regparam,
                )
            end

            result, logdetN = eval_cond(PsiN, Psi, M, S_kGram, S_tilde, ps, zeta, logdet, kappa)
            #println([iter, cnt, zeta, logdetN])

            if result == 1
                break
            end

            if cnt < ninner - 1
                zeta *= c
            elseif cnt == ninner - 1 # safe step size is not actually safe
                min_eig = sum([minimum(eigvals(Psi[k])) for k = 1:K])
                zeta = min_eig ^ 2 / 2
            else # cnt == ninner
                @warn "Backtracking fails. Use more inner iterations / smaller c."
            end
        end

        # Compute convergence
        dPsi = [PsiN[k] - Psi[k] for k = 1:K]
        nrm = norm([dPsi[k] - I * mean(diag(dPsi[k])) for k = 1:K])^2
        out[iter] = [fun(S_kGram, PsiN), reg(PsiN, regcoef), nrm, PsiN, zeta]
        if iter > 3 && all([out[i][3] <= tol for i in iter-3:iter])
            break
            # Stopping criterion based on fun val doesn't work. f keeps decreasing.
            # The problem is the inner product in f doesn't work properly.
            # Is there a reason why Kristjan uses nrm as the criterion?
            # if abs(out[iter][1] - out[iter-1][1]) < tol
            #     break
            # end
        end

        # Set next initial step
        ~, ~, MsN = project_ksum(PsiN, ps)
        zeta = barzilai_borwein(PsiN, Psi, MsN, M, ps)

        # Update iterate
        M .= MsN  # (copy, deepcopy, vs .=) Use === to check if identical
        Psi .= PsiN
        logdet = logdetN

        #println(iter, ": ", [out[iter][1], out[iter][2]])
        #display(heatmaps(PsiN))
        #println("-----------------------------")
    end
    return Psi, out
end

function barzilai_borwein(PsiN, Psi, MsN, M, ps)
    Num = 0
    Den = 0
    K = size(PsiN, 1)
    dgI = zeros(K)
    for k = 1:K
        dgI[k] = mean(diag((PsiN[k] - Psi[k])))
        Num = Num + prod(ps[[1:k-1; k+1:K]]) * sum((PsiN[k] - Psi[k] - I * dgI[k]) .^ 2)
        dK = PsiN[k] - Psi[k]
        for j = 1:K
            dKK = M[j] - MsN[j]
            if k == j
                Den = Den + prod(ps[[1:k-1; k+1:K]]) * sum(sum(dK .* dKK))
            else
                Den = Den + tr(dK) * tr(dKK) * prod(ps) / (ps[k] * ps[j])
            end
        end
    end
    Num = Num + mean(dgI)^2 * prod(ps)
    zeta = Num / Den
    if Den == 0
        zeta = 0
    elseif abs(zeta) > 1e2
        zeta = 1e2 * sign(zeta)
    elseif zeta < 0
        zeta = 1
    end
    return zeta
end

function eval_cond(PsiN, Psi, M, S_kGram, S_tilde, ps, zeta, logdet, kappa)
    K = length(ps)
    p = prod(ps)
    ms = [p / ps[k] for k = 1:K]
    ergNvec = zeros(p, 1) # eigval of the Kronecker sum of PsiN
    for k = 1:K
        ergNvec = ergNvec + repeat(
            kron(eigvals(PsiN[k]), ones(prod(ps[k+1:end]),1)),
            prod(ps[1:k-1]),
            1,
        ) # repmat is eqv to
    end
    #ergNvec = reduce(⊕, [eigvals(PsiN[k]) for k in 1:K])

    if minimum(ergNvec) <= 0
        logdetN = NaN
        result = 0
        return result, logdetN
    end

    logdetN = -sum(log.(ergNvec))
    if maximum(ergNvec) > kappa
        result = 0
        return result, logdetN
    end

    fNew = logdetN + sum([ms[k] * sum(S_kGram[k] .* PsiN[k]) for k = 1:K])
    # Ψ = PsiN[1] ⊕ PsiN[2]
    # SSSS = S_kGram[1] ⊕ S_kGram[2]
    # my_fNew = -log(det(Ψ)) + sum(Ψ .* SSSS)
    # @bp # sum(Ψ .* SSSS) != sum([ms[k] * sum(S_kGram[k] .* PsiN[k]) for k = 1:K])
    # println([fNew, fun(S_kGram, PsiN), my_fNew])

    dPsi = [PsiN[k] - Psi[k] for k = 1:K]
    tau = [mean(diag(dPsi[k])) for k = 1:K] # dgI, tr(D_k)/d_k
    Num = sum([ms[k] * norm(dPsi[k] - tau[k] * I)^2 for k = 1:K])
    Num += mean(tau)^2 * p  # Norm of change
    Q = logdet# + sum([ms[k] * sum(S_kGram[k].*Psi[k]) for k = 1:K]) + Num*K/(2*zeta)
    for k = 1:K #
        Q += sum(S_kGram[k] .* Psi[k]) * ms[k] + Num / (2 * zeta)
        for j = 1:K
            if j == k
                Q += sum((S_tilde[k] - M[k]) .* dPsi[k]) * prod(ps) / ps[k]
            else
                Q -= tr(M[k] - S_tilde[k]) * tr(dPsi[j]) * prod(ps) / (ps[k]*ps[j])
            end
        end
    end

    result = fNew <= Q

    return result, logdetN
end

function fun(S_kGram, Ψs)
    # - logdet(Ω) + < Ω,S >
    ps = [size(S, 1) for S in S_kGram]
    K = length(ps)
    p = prod(ps)
    T = zeros(ps...)
    for (k, Ψ) in enumerate(Ψs)
        dims = ones(Int, length(ps))
        dims[k] = ps[k]
        T .+= reshape(eigvals(Ψ), dims...)
    end
    # 2x faster than this:
    # T = diag(reduce(⊕, [Diagonal(eigvals(Ψs[k])) for k in 1:K]))
    logdet = sum(log.(T))

    trace = sum([sum(S_kGram[k] .* Ψs[k]) * p / ps[k] for k in 1:K])

    return -logdet + trace
end

function reg(Ψs, ρs)
    ps = [size(Ψ, 1) for Ψ in Ψs]
    p = prod(ps)
    result = 0
    for k in 1:length(Ψs)
        mk = p ÷ ps[k]
        offdiag = Ψs[k] - diagm(diag(Ψs[k]))
        result += mk * ρs[k] * sum(abs.(offdiag))
    end
    return result
end


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

@testset "ParaxialLens" begin
    @testset "Virtual point" begin
        lens = ParaxialLensRect(10.0,100.0,100.0,[0.0,0.0,1.0],[0.0,0.0,0.0])
        displaypoint = [0.0,0.0,-8.0]
        r1 = Ray(displaypoint,[1.0,0.0,1.0])
        intsct = OpticSim.surfaceintersection(lens,r1)
        intsctpt = point(intsct)
        # compute the refracted ray and project this backwards to its intersection with the optical axis. This should be the same position as returned by virtualpoint()
        refrac,_,_ = OpticSim.processintersection(OpticSim.interface(lens),OpticSim.point(intsct),OpticSim.normal(lens),OpticalRay(r1,1.0,.55),OpticSim.TEMP_REF,OpticSim.PRESSURE_REF,false)
        # compute intersection of ray with optical axis. this should match the position of the virtual point
        slope = refrac[1]/refrac[3]
        virtptfromslope = [0.0,0.0,-intsctpt[1]/slope]
        # insert test virtptfromslope == point(virtualpoint(lens,displaypoint))
        @test isapprox(virtptfromslope,point(virtualpoint(lens,displaypoint)))
    end
end


text:
######
# This file is part of the MomentArithmetic.jl package, for performing arithmetic operations
# between moments of uncertain numbers, when the moments and the uncertain numbers are only
# partially known.
#
#   Definition of Moment Arithmetic
#
#           University of Liverpool,
#           Institute for Risk and Unertainty
#
#                                           Authors: Ander Gray, Scott Ferson
#                                           Email:  ander.gray@liverpool.ac.uk
######

##
#   Rowe mean
##

function roweFormula(EX, VX, LX , GX, t :: Function)

    P = 1 / (1 + (EX - LX)^2/ VX)
    Q = 1 / (1 + (GX - EX)^2/ VX)

    MIN = P * t(LX) + (1 - P) * t(EX + VX/(EX - LX))
    MAX = Q * t(GX) + (1 - Q) * t(EX + VX/(EX - GX))

    return env(MIN, MAX)

end

function rowe(x :: AbstractMoment, t :: Function; Nsub = Nsub[1])

    if !hasIntMoments(x); return roweNoSub(x,t); end

    EXs = x.mean; VXs = x.var;
    LX = left(x); GX = right(x);

    if typeof(EXs) <: Interval;  EXs = split(EXs, Nsub); end    # Sub-intervalise
    if typeof(VXs) <: Interval;  VXs = split(VXs, Nsub); end

    outs = [roweFormula(EX, VX, LX, GX, t) for EX in EXs, VX in VXs]

    return hull(outs[:])        # Same as env

end

# No sub-intervalisation
function roweNoSub(x :: AbstractMoment, t :: Function)

    EX = x.mean; VX = x.var;
    LX = left(x); GX = right(x);

    return roweFormula(EX, VX, LX, GX, t)

end

###
#   Rowe Variance
##

function rowevarFormula(EX, VX, LX , GX, t :: Function, invT :: Function)

    v = invT(roweFormula(EX, VX, LX , GX, t))

    Lv = left(v); Gv = right(v);

    MIN = ((t(Lv) - t(LX))/(Lv - LX))^2 * (VX + (Lv - EX)^2)
    MAX = ((t(Gv) - t(GX))/(Gv - GX))^2 * (VX + (Gv - EX)^2)

    if isnan(MIN); MIN =0; end
    if isnan(MAX); MAX =0; end

    return env(MIN, MAX)

end

function rowevar(x :: AbstractMoment, t :: Function, invT :: Function; Nsub = Nsub[1])

    if !hasIntMoments(x); return rowevarNoSub(x,t, invT); end

    EXs = x.mean; VXs = x.var;
    LX = left(x); GX = right(x);

    if typeof(EXs) <: Interval;  EXs = split(EXs, Nsub); end    # Sub-intervalise
    if typeof(VXs) <: Interval;  VXs = split(VXs, Nsub); end

    outs = [rowevarFormula(EX, VX, LX, GX, t, invT) for EX in EXs, VX in VXs]

    out = hull(outs[:])
    if isnan(out); out = interval(0, Inf); end
    return out

end

# No sub-intervalisation
function rowevarNoSub(x :: AbstractMoment, t :: Function, invT :: Function)

    EX = x.mean; VX = x.var;
    LX = left(x); GX = right(x);

    return rowevarFormula(EX, VX, LX, GX, t, invT)

end

###
#   Univariate transformations
###
-(x :: AbstractMoment) = Moments(-x.mean, x.var, -x.range)

function reciprocal(x :: AbstractMoment)

    if 0 ∈ x.range
        throw(ArgumentError("Tried to evaluate 1/X, with 0 ∈ X. Interval range is [-∞, ∞], moments not defined"))
    end

    f(x) = 1/x
    invF(x) = 1/x

    MEAN = rowe(x, f)
    VAR = rowevar(x, f, invF)
    RANGE = f(x.range)

    return Moments(MEAN, VAR, RANGE)

end

reciprocal(x :: Real) = 1/x

inv(x :: AbstractMoment) = reciprocal(x)

function exp(x :: AbstractMoment)

    MEAN = rowe(x, exp)
    VAR = rowevar(x, exp, log)
    RANGE = exp(x.range)

    MEAN2 = rowe(2 * x, exp)
    homespun_var = MEAN2 - MEAN^2

    VAR = VAR ∩ homespun_var

    return Moments(MEAN, VAR, RANGE)

end

#^(x :: AbstractMoment) = throw(ArgumentError("rowe exponentiate must be implemented"))

function ^(x :: AbstractMoment, a :: Integer)

    if a == 0;
        return Moments(1, 0, interval(1));
    end

    if a == 1; return x; end

    if a < 0;
        if 0 ∈ x.range
            throw(ArgumentError("Tried to evaluate 1/0, when 0 ∈ X when performing X^a where a < 0. Interval range is [-∞, ∞], Moments not defined"))
        end
        return reciprocal(x)^a;
    end

    if isodd(a) && 0 > x.range.lo && 0 < x.range.hi;
        return multFrechet(x^(a-1), x);
    end

    f(x) = x^a
    invF(x) = x^(1/a)

    MEAN = rowe(x, f)
    VAR = rowevar(x, f, invF)
    RANGE = f(x.range)

    if a == 2
        var_hs = homespun_var_square(x)
        VAR = VAR ∩ var_hs
    end

    return Moments(MEAN, VAR, RANGE)
end

function homespun_var_square( x :: Moments )
    that = x^4
    this = that.mean - (x.mean^2 + x.var)^2
    return max(0,this)
end

function ^(x :: AbstractMoment, a :: Float64)

    if isinteger(a); return x^Integer(a); end

    if x.range.lo < 0;
        throw(ArgumentError("Cannot take the power of a negative number to a real number. Yields imaginary numbers, for which moments are currently undefined"))
    end

    f(x) = x^a
    invF(x) = x^(1/a)

    MEAN = rowe(x, f)
    VAR = rowevar(x, f, invF)
    RANGE = f(x.range)

    return Moments(MEAN, VAR, RANGE)
end

function ^(x :: AbstractMoment, a :: Interval)

    if isinteger(a); return x^Integer(a); end

    if x.range.lo < 0;
        throw(ArgumentError("Cannot take the power of a negative number to a real number. Yields imaginary numbers, for which moments are currently undefined"))
    end

    x1 = x^(a.lo)
    x2 = x^(a.hi)

    out = env(x1, x2)

    unitMom = Moments(1,0,interval(1))

    if 0 ∈ a; out = env(out, unitMom); end
    if 1 ∈ x.range; out = env(out, unitMom); end

    return out
end

function sqrt(x :: AbstractMoment)

    f(x) = sqrt(x)
    invF(x) = x^2

    MEAN = rowe(x, f)
    VAR = rowevar(x, f, invF)
    RANGE = f(x.range)

    homespun_var = x.mean - MEAN^2
    VAR = VAR ∩ homespun_var;

    return Moments(MEAN, VAR, RANGE)
end

function ln(x :: AbstractMoment)

    f(x) = log(x)
    invF(x) = exp(x)

    MEAN = rowe(x, f)
    VAR = rowevar(x, f, invF)
    RANGE = f(x.range)

    return Moments(MEAN, VAR, RANGE)
end

function log(x = missing :: AbstractMoment, base = ℯ :: Real)
    if base == ℯ; return ln(x); end

    f(x) = log(base, x)
    invF(x) = base^x

    MEAN = rowe(x, f)
    VAR = rowevar(x, f, invF)
    RANGE = f(x.range)

    return Moments(MEAN, VAR, RANGE)
end

function abs(x :: AbstractMoment)
    if x.range >= 0;
        return x
    elseif x.range < 0;
        EY = -x.mean;
    else
        sqvar = sqrt(x.var)
        EY1 = abs(x.mean)
        EY2 = EY1 + sqvar * (pi - atan(EY1/sqvar));
        EY = env(EY1, EY2)
    end

    Yvar = max(0, x.mean^2 + x.var - EY);
    Yrange = abs(x.range)

    return Moments(EY, Yvar, Yrange)

end

###
#   Scalar arithmetic
###
+(x :: AbstractMoment, y :: Real) = Moments(x.mean + y, x.var, x.range + y)
+(x :: Real, y:: AbstractMoment) = y + x;
-(x :: AbstractMoment, y :: Real) = x + (-y);
-(x :: Real, y :: AbstractMoment) = x + (-y);

*(x :: AbstractMoment, y :: Real) = Moments(x.mean * y, x.var * y^2, x.range * y)
*(x :: Real, y:: AbstractMoment) = y * x;
/(x :: AbstractMoment, y :: Real) = x * (1/y);
/(x :: Real, y :: AbstractMoment) = x * reciprocal(y);

###
#   Moment Arithmetic
###

###
# Independent
###
function sumIndep(x :: AbstractMoment, y :: AbstractMoment)

    zMean  = x.mean  + y.mean;
    zVar   = x.var   + y.var;
    zRange = x.range + y.range;

    # Add dependent variables to dict
    return Moments(zMean, zVar, zRange);
end

function subIndep(x :: AbstractMoment, y :: AbstractMoment)

    zMean  = x.mean  - y.mean;
    zVar   = x.var   + y.var;
    zRange = x.range - y.range;

    # Add dependent variables to dict
    return Moments(zMean, zVar, zRange);
end

function multIndep(x :: AbstractMoment, y :: AbstractMoment)

    zMean  = x.mean  * y.mean;
    zVar   = (x.mean^2 * y.var) + (y.mean^2 * x.var) + (x.var * y.var);     # No subint required
    zRange = x.range * y.range;

    # Add dependent variables to dict
    return Moments(zMean, zVar, zRange);

end

function divIndep(x :: AbstractMoment, y :: AbstractMoment)
    multIndep(x,1/y);
end

function minIndep(x :: AbstractMoment, y :: AbstractMoment)

    if x.range.hi < y.range.lo
        zMean = x.mean;
        zVar = x.var;
        zRange = min(x.range, y.range)
        return Moments(zMean, zVar, zRange)
    end

    if y.range.hi < x.range.lo
        zMean = y.mean;
        zVar = y.var;
        zRange = min(x.range, y.range)
        return Moments(zMean, zVar, zRange)
    end

    return minFrechet(x, y)
end

function maxIndep(x :: AbstractMoment, y :: AbstractMoment)

    if x.range.hi < y.range.lo
        zMean = y.mean;
        zVar = y.var;
        zRange = max(x.range, y.range)
        return Moments(zMean, zVar, zRange)
    end

    if y.range.hi < x.range.lo
        zMean = x.mean;
        zVar = x.var;
        zRange = max(x.range, y.range)
        return Moments(zMean, zVar, zRange)
    end

    return maxFrechet(x, y)
end

+(x :: AbstractMoment, y :: AbstractMoment) = sumFrechet(x,y);
-(x :: AbstractMoment, y :: AbstractMoment) = subFrechet(x,y);
*(x :: AbstractMoment, y :: AbstractMoment) = multFrechet(x,y);
/(x :: AbstractMoment, y :: AbstractMoment) = divFrechet(x,y);
min(x :: AbstractMoment, y :: AbstractMoment) = minFrechet(x,y);
max(x :: AbstractMoment, y :: AbstractMoment) = maxFrechet(x,y);

##
# With Intervals (Moments with just ranges)
##
+(x :: AbstractMoment, y :: Interval) = sumFrechet(x, Moments(missing, missing, y));
-(x :: AbstractMoment, y :: Interval) = subFrechet(x, Moments(missing, missing, y));
*(x :: AbstractMoment, y :: Interval) = multFrechet(x, Moments(missing, missing, y));
/(x :: AbstractMoment, y :: Interval) = divFrechet(x,  Moments(missing, missing, y));
min(x :: AbstractMoment, y :: Interval) = minFrechet(x,Moments(missing, missing, y));
max(x :: AbstractMoment, y :: Interval) = maxFrechet(x,Moments(missing, missing, y));

+(x :: Interval, y :: AbstractMoment) = sumFrechet(Moments(missing, missing, x), y);
-(x :: Interval, y :: AbstractMoment) = subFrechet(Moments(missing, missing, x), y);
*(x :: Interval, y :: AbstractMoment) = multFrechet(Moments(missing, missing, x),y);
/(x :: Interval, y :: AbstractMoment) = divFrechet(Moments(missing, missing, x), y);
min(x :: Interval, y :: AbstractMoment) = minFrechet(Moments(missing, missing, x),y);
max(x :: Interval, y :: AbstractMoment) = maxFrechet(Moments(missing, missing, x),y);

###
#   Dependence unknown
###
function sumFrechet(x :: AbstractMoment, y :: AbstractMoment)

    zMean = x.mean + y.mean

    zVarLb = (sqrt(x.var) - sqrt(y.var))^2
    zVarUb = (sqrt(x.var) + sqrt(y.var))^2

    zVarLb = max(zVarLb, 0)

    zVar = env(zVarLb, zVarUb);

    zRange = x.range + y.range;

    return Moments(zMean, zVar, zRange);
end

function subFrechet(x :: AbstractMoment, y :: AbstractMoment)

    zMean = x.mean - y.mean

    zVarLb = (sqrt(x.var) - sqrt(y.var))^2
    zVarUb = (sqrt(x.var) + sqrt(y.var))^2

    zVarLb = max(zVarLb, 0)

    zVar = env(zVarLb, zVarUb);

    zRange = x.range - y.range;

    return Moments(zMean, zVar, zRange);
end

###
#   Frechet multiplication. Not sub-intervalising required.
###
function multFrechet(x :: AbstractMoment, y :: AbstractMoment)

    EX = x.mean; EY = y.mean;
    VX = x.var; VY = y.var;

    zMeanLb = EX * EY - sqrt(VX * VY);
    zMeanUb = EX * EY + sqrt(VX * VY);
    zMean = env(zMeanLb, zMeanUb);

    x2 = x^2; y2 = y^2

    E11 = cov(x2,y2) + (x2.mean * y2.mean);
    E22 = (cov(x,y) + (x.mean * y.mean)) ^2

    zVar = E11 - E22

    zRange = x.range * y.range;

    Moments(zMean, zVar, zRange)

end

# Requires alot of subintervaling for Gvar to have an effect
function multFrechetOld(x :: AbstractMoment, y :: AbstractMoment; Nsub = Nsub[1])

    EX = x.mean; EY = y.mean;
    VX = x.var; VY = y.var;

    zMeanLb = EX * EY - sqrt(VX * VY);
    zMeanUb = EX * EY + sqrt(VX * VY);
    zMean = env(zMeanLb, zMeanUb);

    if typeof(EX) <: Interval && !isvacuous(EX);  EX = split(EX, Nsub); end    # Sub-intervalise
    if typeof(EY) <: Interval && !isvacuous(EY);  EY = split(EY, Nsub); end

    zVar = [Gvar(Moments(ex ,VX ,x.range), Moments(ey, VY, y.range)) for ex in EX, ey in EY]

    zVar = hull(zVar[:])

    zRange = x.range * y.range;

    return Moments(zMean, zVar, zRange);
end

divFrechet(x :: AbstractMoment, y :: AbstractMoment) = multFrechet(x, 1/y)  #Not best possible

function minFrechet(x :: AbstractMoment, y :: AbstractMoment)

    zRange = min(x.range, y.range)
    zVar = env(0, max(x.var, y.var))

    EX = x.mean; EY = y.mean;

    LX = left(x.range); GX = right(x.range)
    LY = left(y.range); GY = right(y.range)

    zMean = Bmin(EX, LX, GX, EY, LY, GY);
    return Moments(zMean, zVar, zRange)

end

function maxFrechet(x :: AbstractMoment, y :: AbstractMoment)

    zRange = max(x.range, y.range)
    zVar = env(0, max(x.var, y.var))

    EX = x.mean; EY = y.mean;

    LX = left(x.range); GX = right(x.range)
    LY = left(y.range); GY = right(y.range)

    zMean = Bmax(EX, LX, GX, EY, LY, GY);
    return Moments(zMean, zVar, zRange)

end

###
#   Perfect Arithmetic
###

function sumPerfect(x :: Moments, y :: Moments)

    meanZ = x.mean + y.mean
    varZ  = (sqrt(x.var) + sqrt(y.var))^2
    rangeZ = x.range + y.range
    return Moments(meanZ, varZ, rangeZ)
end

function subPerfect(x :: Moments, y :: Moments)

    meanZ = x.mean - y.mean
    varZ  = (sqrt(x.var) - sqrt(y.var))^2
    rangeZ = x.range - y.range
    return Moments(meanZ, varZ, rangeZ)
end

multPerfect(x, y) = multCor(x, y, 1)

###
#   Opposite Arithmetic
###

function sumOpposite(x :: Moments, y :: Moments)

    meanZ = x.mean + y.mean
    varZ  = (sqrt(x.var) - sqrt(y.var))^2
    rangeZ = x.range + y.range
    return Moments(meanZ, varZ, rangeZ)
end

function subOpposite(x :: Moments, y :: Moments)

    meanZ = x.mean - y.mean
    varZ  = (sqrt(x.var) + sqrt(y.var))^2
    rangeZ = x.range - y.range
    return Moments(meanZ, varZ, rangeZ)
end

###
#   known Dependence, or interval Dependence
###

function sumCor(x :: Moments, y :: Moments, corr :: Interval, nSub = 10)

    xVar = x.var; yVar = y.var; R = corr;

    if typeof(x.var) <: Interval; xVar = mince(xVar, nSub); end
    if typeof(y.var) <: Interval; yVar = mince(yVar, nSub); end
    #if typeof(corr) <: Interval; R = mince(R, nSub); end

    varZ = [xs + ys + 2 * sqrt(xs) * sqrt(ys) * corrs for xs in xVar, ys in yVar, corrs in R]

    meanZ = x.mean + y.mean
    rangeZ = x.range + y.range

    varZ = hull(varZ[:]);

    return Moments(meanZ, varZ, rangeZ)
end

sumCor(x :: Moments, y :: Moments, corr :: Real, nSub = 10) = sumCor(x,y, interval(corr), nSub)

function subCor(x :: Moments, y :: Moments, corr :: Interval, nSub = 10)

    xVar = x.var; yVar = y.var; R = corr;

    if typeof(x.var) <: Interval; xVar = mince(xVar, nSub); end
    if typeof(y.var) <: Interval; yVar = mince(yVar, nSub); end
    #if typeof(corr) <: Interval; R = mince(R, nSub); end

    varZ = [xs + ys - 2 * sqrt(xs) * sqrt(ys) * corrs for xs in xVar, ys in yVar, corrs in R]

    meanZ = x.mean - y.mean
    rangeZ = x.range - y.range

    varZ = hull(varZ[:]);

    return Moments(meanZ, varZ, rangeZ)
end

subCor(x :: Moments, y :: Moments, corr :: Real, nSub = 10) = subCor(x,y, interval(corr), nSub)

function multCor(x :: Moments, y :: Moments, corr :: Interval)

    EX = x.mean; EY = y.mean;
    VX = x.var; VY = y.var;

    zMean = EX * EY + corr * sqrt(VX * VY);

    x2 = x^2; y2 = y^2

    E11 = cov(x2,y2) + (x2.mean * y2.mean);
    E22 = (corr * sqrt(VX * VY) + (x.mean * y.mean)) ^2

    zVar = E11 - E22

    zRange = x.range * y.range;

    return Moments(zMean, zVar, zRange)
end

multCor(x :: Moments, y :: Moments, corr :: Real) = multCor(x,y, interval(corr))

function divCor(x :: Moments, y :: Moments, corr :: Interval)
    return multCor(x, 1/y, -1 * corr)
end

divCor(x :: Moments, y :: Moments, corr :: Real) = divCor(x,y, interval(corr))

###
#   Using Covariances
###

function sumCov(x :: Moments, y :: Moments, cov :: Interval)

    xVar = x.var; yVar = y.var;

    varZ = xVar + yVar + 2 * cov;

    meanZ = x.mean + y.mean
    rangeZ = x.range + y.range

    return Moments(meanZ, varZ, rangeZ)
end

sumCov(x :: Moments, y :: Moments, cov :: Real) = sumCov(x,y, interval(cov))

function subCov(x :: Moments, y :: Moments, cov :: Interval)

    xVar = x.var; yVar = y.var;

    meanZ = x.mean - y.mean

    varZ = xVar + yVar - 2 * cov;

    rangeZ = x.range - y.range

    return Moments(meanZ, varZ, rangeZ)
end

subCov(x :: Moments, y :: Moments, cov :: Real) = subCov(x,y, interval(cov))

function multCov(x :: Moments, y :: Moments, Cov :: Interval)

    EX = x.mean; EY = y.mean;

    zMean = EX * EY + Cov;

    x2 = x^2; y2 = y^2

    E11 = cov(x2,y2) + (x2.mean * y2.mean);
    E22 = (Cov + (x.mean * y.mean))^2

    zVar = E11 - E22

    zRange = x.range * y.range;

    return Moments(zMean, zVar, zRange)
end

multCov(x :: Moments, y :: Moments, cov :: Real) = multCov(x,y, interval(cov))

###
#   Bertsimas mean for min and max. Do not require sub-intervalisation due to env and intersect
###

function Bmin(EX, LX, GX, EY, LY, GY)

    B1 = env(min(EX, EY), min(LX, LY))
    B2 = EX + EY - env(max(EX, EY), max(GX, GY))

    return intersect(B1,B2)
end

function Bmax(EX, LX, GX, EY, LY, GY)

    B1 = env(max(EX, EY), max(GX, GY))
    B2 = EX + EY - env(min(EX, EY), min(LX, LY))

    return intersect(B1,B2)
end

###
#   Goodman formula for variance under unknown dependence
###

#   Looks like a complidated repeat variables problem, but only EX and EY need sub-ints
function Gvar(x :: AbstractMoment, y :: AbstractMoment)

    EX = x.mean; VX = y.var
    EY = y.mean; VY = y.var

    EXY = multFrechetMean(x, y)
    EX2Y = multFrechetMean(x^2, y)
    EXY2 = multFrechetMean(x, y^2)
    EX2Y2 = multFrechetMean(x^2, y^2)

    EX2 = (x^2).mean
    EY2 = (y^2).mean

    E11 = EXY - EX*EY
    E12 = EX2Y - EX2*EY + 2*EX^2*EY - 2*EX*EXY
    E21 = EXY2 - EX*EY2 + 2*EX*EY^2 - 2*EY*EXY
    E22 = -3*EX^2*EY^2 + EX2*EY^2 + EX^2*EY2 + 4*EX*EY*EXY - 2*EY*EX2Y - 2*EX*EXY2 + EX2Y2

    VXY = EX^2*VX + EY^2*VY + 2*EX*EY*E11 + 2*EX*E12 + 2*EY*E21 + E22 - E11^2

    #VXY2 = EX^2*VX + EY^2*VY + 2*EX*EX2Y - 2*EY*EX2Y + 4*EX*EY^3 - 6*EX^2*EY^2 + 4*EX^3*EY - EXY^2 + 8*EX*EY*EXY - 4*EX^2*EXY - 4*EY^2*EXY -2*EX*EY*EX2 + EY^2*EX2 - 2*EX*EY*EY2 + EX^2*EY2 + 2*EY*EXY2 - 2*EX*EXY2 + EX2Y2

    return max(VXY, 0)

end

#=

    E11 = EXY - EX*EY
    E12 = EX2Y - EX2*EY + 2*EX^2*EY - 2*EX*EXY
    E21 = EXY2 - EX*EY2 + 2*EX*EY^2 - 2*EY*EXY
    E22 = -3*EX^2*EY^2 + EX2*EY^2 + EX^2*EY2 + 4*EX*EY*EXY - 2*EY*EX2Y - 2*EX*EXY2 + EX2Y2

    VXY = EX^2*VX + EY^2*VY + 2*EX*EY*E11 + 2*EX*E12 + 2*EY*E21 + E22 - E11^2

    A = EX; B = EY;
    C = EXY; D = EX2;
    E = EY2; F = EX2Y;
    G = EXY2; H = EX2Y2

    2 * A * B * (C - A*B) + 2 * A * (F -D*B +2 * A^2 * B - 2 * A * C) + 2 * B * (G - A * E + 2 * A * B^2 - 2 * B * C) + (-3 * A^2 * B^2 + D * B^2 + A^2 * E + 4*A*B*C -2*B*F - 2*A*G + H) - (C - A*B)^2

    simplify gives:

    2AF-2BF+4AB^3-6A^2B^2+4A^3B-C^2+8ABC-4A^2C-4B^2C-2ABD+B^2D-2ABE+A^2E+2BG-2AG+H

    2*EX*EX2Y - 2*EY*EX2Y + 4*EX*EY^3 - 6*EX^2*EY^2 + 4*EX^3*EY - EXY^2 + 8*EX*EY*EXY - 4*EX^2*EXY - 4*EY^2*EXY -2*EX*EY*EX2 + EY^2*EX2 - 2*EX*EY*EY2 + EX^2*EY2 + 2*EY*EXY2 - 2*EX*EXY2 + EX2Y2
=#

###
#   Woodpile
###

#=
###
#   Bertsimas with subintervalistion (does not require)
###
function minFrechetSub(x :: Moments, y ::Moments)

    zRange = min(x.range, y.range)
    zVar = env(0, max(x.var, y.var))

    EX = x.mean; EY = y.mean;

    LX = left(x.range); GX = right(x.range)
    LY = left(y.range); GY = right(y.range)

    if !hasIntMean(x) && !hasIntMean(y);
        zMean = Bmin(EX, LX, GX, EY, LY, GY);
        return Moments(zMean, zVar, zRange)
    end

    if typeof(EX) <: Interval;  EX = Split(EX, Nsub); end
    if typeof(EY) <: Interval;  EY = Split(EY, Nsub); end

    zMeans = [Bmin(EX, LX, GX, EY, LY, GY) for EX in EX, EY in EY];

    zMean = hull(zMeans[:])

    return Moments(zMean, zVar, zRange)

end

F(xmean, ymean, xvar, yvar) = (xmean^2 * yvar) + (ymean^2 * xvar) + (xvar * yvar)

xmean = interval(0,100)
xvar = interval(0,3)
ymean = interval(0,100)
yvar = interval(0,3)

xs = Split(xvar, 10)
ys = Split(yvar, 10)

outs  = [F(xmean, ymean, xvar, yvar) for xvar in xs, yvar in ys]
hull(outs[:])

function multFrechet(x :: AbstractMoment, y :: AbstractMoment; sub = false, Nsub = Nsub[1])

    EX = x.mean; EY = y.mean;
    VX = x.var; VY = y.var;

    zMeanLb = EX * EY - sqrt(VX * VY);
    zMeanUb = EX * EY + sqrt(VX * VY);
    zMean = env(zMeanLb, zMeanUb);

    if sub
        if typeof(VX) <: Interval && !isvacuous(VX);  VX = split(VX, Nsub); end    # Sub-intervalise
        if typeof(VY) <: Interval && !isvacuous(VY);  VY = split(VY, Nsub); end
    end

    COR = interval(-1, 1)           # cor = [-1, 1] for Frechet

    zVar = [(VX + EX^2) * (VY + EY^2) + COR*VX*VY - (COR * sqrt(VX) * sqrt(VY) + EY * EX)^2 for VX in VX, VY in VY]
    zVar = hull(zVar[:])

    zRange = x.range * y.range;

    Moments(zMean, zVar, zRange)

end

=#


text:
"""
Collection of logic for how response data is parsed into concrete data types.
"""
module Parsers

using Stonks.Models: AbstractStonksRecord

export AbstractContentParser, CSVParser, JSONParser, parse_content

"""
Abstract type to be subclassed by any type of parser.
"""
abstract type AbstractContentParser end

"""
A wrapper type for a function implementing `parse_content` interface for JSON content.
"""
struct JSONParser{F} <: AbstractContentParser
  func::F
end

"""
A wrapper type for a function implementing `parse_content` interface for CSV content.
"""
struct CSVParser{F} <: AbstractContentParser
  func::F
end

function parse_content end

"""
    parse_content(p::AbstractContentParser, c::AbstractString; kwargs...) -> Vector{T<:AbstractStonksRecord}

Interface for parsing raw data into a data model.
All main functions responsible for getting data from APIs will call `parse_content` with a concrete subtype of `AbstractContentParser`.

### Example - Create a custom type and parser
```julia
using Dates, JSON3, Stonks

# Assume we have the following raw data
julia> content
{
  "name": "Inflation - US Consumer Prices",
  "interval": "annual",
  "unit": "percent",
  "data": [
    {"date": "2020-01-01", "value": "1.2335"},
    {"date": "2019-01-01", "value": "1.8122"}
  ]
}

# Create your custom type
struct MacroIndicator <: AbstractStonksRecord
  name::String
  date::Date
  value::Number
end

# Define a function with the same signature as parse_content.
function parse_inflation_data(content::AbstractString; kwargs...)
  js = JSON3.read(content)
  return map(x ->
    MacroIndicator("Inflation", tryparse(Date, x[:date]), tryparse(Float64, x[:value])),
  js["data"])
end

# Wrap it around a concrete subtype of AbstractContentParser
my_parser = Stonks.Parsers.JSONParser(parse_inflation_data)
julia> parse_content(my_parser, content)
2-element Vector{MacroIndicator}:
 MacroIndicator("Inflation", Date("2020-01-01"), 1.2335)
 MacroIndicator("Inflation", Date("2019-01-01"), 1.8122)
```
"""
parse_content(p::AbstractContentParser, c::AbstractString; kwargs...) = p.func(c; kwargs...)

include("parsers/utils.jl")
include("parsers/parser_yahoo.jl")
include("parsers/parser_alphavantage.jl")

YahooPriceParser = JSONParser(parse_yahoo_price)
YahooInfoParser = JSONParser(parse_yahoo_info)
YahooExchangeRateParser = JSONParser(parse_yahoo_exchange_rate)
YahooIncomeStatementParser = JSONParser(parse_yahoo_income_statement)
YahooBalanceSheetParser = JSONParser(parse_yahoo_balance_sheet)
YahooCashflowParser = JSONParser(parse_yahoo_cashflow_statement)
YahooEarningsParser = JSONParser(parse_yahoo_earnings)

AlphavantagePriceParser = JSONParser(parse_alphavantage_price)
AlphavantageInfoParser = JSONParser(parse_alphavantage_info)
AlphavantageExchangeRateParser = JSONParser(parse_alphavantage_exchange_rate)
AlphavantageIncomeStatementParser = JSONParser(parse_alphavantage_income_statement)
AlphavantageBalanceSheetParser = JSONParser(parse_alphavantage_balance_sheet)
AlphavantageCashflowParser = JSONParser(parse_alphavantage_cashflow_statement)
AlphavantageEarningsParser = JSONParser(parse_alphavantage_earnings)

end


text:
###########
# ACTIONS #
###########

function createRegexNode(rule, value, first, last, children)
  return Node("REGEX", children[1].value, first, last, children[1].children, children[1].ruleType)
end
function createActionNode(rule, value, first, last, children)
  return Node("ACTION", children[1].value, first, last, children[1].children, children[1].ruleType)
end
function createRefNode(rule, value, first, last, children)
  return Node("REF", children[1].value, first, last, children[1].children, children[1].ruleType)
end
function createTermNode(rule, value, first, last, children)
  disescapedcontent = replace(children[1].value, "''" => "'")
  return Node("TERM", disescapedcontent, first, last, [], Terminal)
end

###########################################################
# HANDCONSTRUCTION OF THE GRAMMAR TO PARSE GRAMMARSTRINGS #
###########################################################

# legibility
sup = SuppressRule
ref = ReferencedRule
and = AndRule
or  = OrRule

"""
The grammar to parse grammars. For obvious reasons this has to be constructed "by hand" from constructors, as it is needed itself for the parsing of grammar strings. `grammargrammar_string` is a consistency check and parses to grammargrammar - but only when grammargrammar already exists for its parsing.
"""
const grammargrammar = Grammar(Dict{Symbol,Any}(
:start     => and([ sup(ZeroOrMoreRule(ref(:emptyline))), ZeroOrMoreRule("ALLRULES",ref(:ruleline)), sup(OptionalRule(ref(:space))) ],liftchild),

:emptyline => and( sup(ref(:space)), ref(:endofline) ),
:ruleline  => and([ sup(ref(:space)), ref(:rule), ZeroOrMoreRule(sup(ref(:emptyline))) ],liftchild),
:rule      => and("RULE",[ ref(:symbol), sup(ref(:space)), sup(Terminal("=>")), sup(ref(:space)), ref(:definition)]),

:definition=> or([ ref(:double), ref(:single) ]), # from left to right: ' can contain ' => parse order

:single    => and("SINGLE", or([ref(:parenrule),ref(:zeromorerule),ref(:onemorerule),ref(:optionalrule),ref(:notrule),ref(:suppressrule),ref(:regexrule),ref(:term),ref(:refrule)]), OptionalRule(and([sup(ref(:space)),ref(:action)],liftchild),liftchild) ), # only single token rules can have associated actions (-> unique interpretation)
:double    => or([ref(:orrule),ref(:andrule)]), # 'and' groups stronger than 'or'

:parenrule => and("PAREN",[ sup(Terminal('(')), sup(ref(:space)), ref(:definition), sup(ref(:space)), sup(Terminal(')')) ]),
:orrule    => and("OR",[ or(ref(:andrule),ref(:single)), OneOrMoreRule("more",and([ sup(ref(:space)), sup(Terminal('|')), sup(ref(:space)), or(ref(:andrule),ref(:single)) ],liftchild)) ]),
:andrule   => and("AND",[ ref(:single), OneOrMoreRule("more",and([ sup(ref(:space)), sup(Terminal('&')), sup(ref(:space)), ref(:single) ],liftchild)) ]),
:zeromorerule => and("*",[ sup(Terminal("*(")), sup(ref(:space)), ref(:definition), sup(ref(:space)), sup(Terminal(')')) ]),
:onemorerule  => and("+",[ sup(Terminal("+(")), sup(ref(:space)), ref(:definition), sup(ref(:space)), sup(Terminal(')')) ]),
:optionalrule => and("?",[ sup(Terminal("?(")), sup(ref(:space)), ref(:definition), sup(ref(:space)), sup(Terminal(')')) ]),
:notrule      => and("!",[ sup(Terminal("!(")), sup(ref(:space)), ref(:definition), sup(ref(:space)), sup(Terminal(')')) ]),
:suppressrule => and("-",[ sup(Terminal("-(")), sup(ref(:space)), ref(:definition), sup(ref(:space)), sup(Terminal(')')) ]),
:refrule   => ref(:symbol,createRefNode),
:term      => and([ sup(Terminal('\'')), RegexRule(r"([^']|'')+"), sup(Terminal('\'')) ], createTermNode),
:regexrule => and([ sup(Terminal("r(")), RegexRule(r".*?(?=\)(?=r))"), sup(Terminal(")r")) ],createRegexNode), # r(...)r to prevent escaping issues with '(', ')' within the regex, '"' has the same issue but is even weirder because we want to parse a string "... r"..." ..." would then actually have to be input as "... r\"...\" ..." and within the regex "... r\" \\" \" ..."

:action    => and([sup(Terminal('{')), RegexRule(r"[^}]*"), sup(Terminal('}'))],createActionNode),

:space     => RegexRule(r"[ \t]*"),
:endofline => or([Terminal("\r\n"), Terminal('\r'), Terminal('\n'), Terminal(';')]),
:symbol    => RegexRule("SYM",r"[a-zA-Z_][a-zA-Z0-9_]*"),
))

######################################
# TRANSFORMS FOR THIS GRAMMARGRAMMAR #
######################################

togrammar(node, children, ::MatchRule{:default}) = (warn("Unmatched transform: $node"); children)
togrammar(node, children, ::MatchRule{:PAREN}) = children[1]
togrammar(node, children, ::MatchRule{:TERM}) = Terminal(node.value)
togrammar(node, children, ::MatchRule{:REGEX}) = RegexRule(Regex(node.value))
togrammar(node, children, ::MatchRule{:REF}) = ReferencedRule(Symbol(node.value))
togrammar(node, children, ::MatchRule{:ACTION}) = node.value
function togrammar(node, children, ::MatchRule{:SINGLE})
  node = children[1]
  if length(children)==2 # action specification exists
    action = Meta.parse(children[2])
    if isa(action,AbstractString)
      node.name = action
    elseif isa(action,Function)
      node.action = action
    elseif isa(action,Expr)
      node.action = eval(action)
    elseif isa(action,Symbol)
      node.action = getfield(StringParserPEG,action)
    else
      error("Unexpected action type $(typeof(action)): $action")
    end
  end
  return node
end
togrammar(node, children, ::MatchRule{:*}) = ZeroOrMoreRule(children[1])
togrammar(node, children, ::MatchRule{:?}) = OptionalRule(children[1])
togrammar(node, children, ::MatchRule{:+}) = OneOrMoreRule(children[1])
togrammar(node, children, ::MatchRule{:!}) = NotRule(children[1])
togrammar(node, children, ::MatchRule{:-}) = SuppressRule(children[1])
togrammar(node, children, ::MatchRule{:more}) = Vector{Rule}(children)
function togrammar(node, children, ::MatchRule{:AND})
  list = Vector{Rule}()
  push!(list,children[1])
  append!(list,children[2])
  AndRule(list)
end
function togrammar(node, children, ::MatchRule{:OR})
  list = Vector{Rule}()
  push!(list,children[1])
  append!(list,children[2])
  OrRule(list)
end
togrammar(node, children, ::MatchRule{:SYM}) = Symbol(node.value)
togrammar(node, children, ::MatchRule{:RULE}) = (children[1],children[2])
togrammar(node, children, ::MatchRule{:ALLRULES}) = Grammar(Dict(children))

###########
# TESTING #
###########

"""
`grammargrammar` parses `grammargrammar_string` such that when transformed with `transform(togrammar,ast)` again `grammmargrammmar` results. This is a consistency check and allows to understand what happens in the `grammargrammar` construction in a more legible form.
"""
const grammargrammar_string = """
start      => (-(*(emptyline)) & *(ruleline) {"ALLRULES"} & -(?(space))) {liftchild}

emptyline  => -(space) & endofline
ruleline   => ( -(space) & rule & *(-(emptyline)) ) {liftchild}
rule       => ( symbol & -(space) & -('=>') & -(space) & definition ) {"RULE"}

definition => double | single

single     => ( (parenrule | zeromorerule | onemorerule | optionalrule | notrule | suppressrule | regexrule | term | refrule) & ?((-(space) & action) {liftchild}) {liftchild} ) {"SINGLE"}
double     => orrule | andrule

parenrule  => ( -('(') & -(space) & definition & -(space) & -(')') ) {"PAREN"}
orrule     => ( (andrule | single) & +( (-(space) & -('|') & -(space) & (andrule|single)){liftchild} ){"more"} ) {"OR"}
andrule    => (            single  & +( (-(space) & -('&') & -(space) &          single ){liftchild} ){"more"} ) {"AND"}
zeromorerule   => ( -('*(') & -(space) & definition & -(space) & -(')') ){"*"}
onemorerule    => ( -('+(') & -(space) & definition & -(space) & -(')') ){"+"}
optionalrule   => ( -('?(') & -(space) & definition & -(space) & -(')') ){"?"}
notrule        => ( -('!(') & -(space) & definition & -(space) & -(')') ){"!"}
suppressrule   => ( -('-(') & -(space) & definition & -(space) & -(')') ){"-"}
refrule    => symbol {createRefNode}
term       => ( -('''') & r(([^']|'')+)r & -('''') ) {createTermNode}
regexrule  => ( -('r(') & r(.*?(?=\\)(?=r)))r & -(')r') ) {createRegexNode}

action     => (-('{') & r([^}]*)r & -('}')) {createActionNode}

space      => r([ \\t]*)r
endofline  => '\r\n' | '\r' | '\n' | ';'
symbol     => r([a-zA-Z_][a-zA-Z0-9_]*)r {"SYM"}
"""