text:
using Base.Test

importall BioSymbols
importall BioSequences
importall BioAlignments
importall IntervalTrees

using DataStructures
using BufferedStreams
using BioSequences
using BioSymbols
using BioAlignments
using FMIndexes
using IntArrays
using IntervalTrees
using Libz
using Distributions
using Requests

include("../src/types.jl")
include("../src/timer.jl")
include("../src/sgkmer.jl")
include("../src/fmindex_patch.jl")
include("../src/bam.jl")
include("../src/refset.jl")
include("../src/graph.jl")
include("../src/bias.jl")
include("../src/edges.jl")
include("../src/index.jl")
include("../src/record.jl")
include("../src/align.jl")
include("../src/quant.jl")
include("../src/reads.jl")
include("../src/ebi.jl")
include("../src/paired.jl")
include("../src/events.jl")
include("../src/io.jl")
include("../src/diff.jl")

@testset "SG Kmers & Seq" begin
   @test sgkmer(dna"ATG") == BioSequences.DNAKmer(dna"ATG")
   @test sgkmer("ATG") == BioSequences.DNAKmer(dna"ATG")
   @test isa(sgkmer(dna"ATG"), SGKmer{3})
   @test kmer_index( sgkmer(dna"ATG") ) == kmer_index( BioSequences.DNAKmer(dna"ATG") )
   @test Int(kmer_index_trailing( dna"ATG" ))+1 == kmer_index( BioSequences.DNAKmer(dna"ATG") )
   @test Int(kmer_index_trailing( dna"ATGR" )) == 0
   seq = dna""
   for i in 1:32
      seq = i % 2 == 0 ? seq * dna"C" : seq * dna"T"
      curkmer = sgkmer(seq)
      @test isa( curkmer, SGKmer{i} )
   end
end

@testset "Bias Models" begin
   # default 'empty' model
   def = zero(DefaultCounter)
   @test get(def) == 0.0
   def = one(DefaultCounter)
   @test get(def) == 1.0
   push!( def, 1.0 )
   @test get(def) == 2.0
   gc_adjust!( def, DefaultBiasMod() )
   primer_adjust!( def, DefaultBiasMod() )

   function randdna(n)
      return BioSequence{DNAAlphabet{2}}(rand([DNA_A, DNA_C, DNA_G, DNA_T], n))
   end

   function test_uniform_bias!( mod::B ) where B <: BiasModel
      for i in 1:5000000
         primer_count!( mod, randdna(36) )
      end
      primer_normalize!( mod )
      for i in 1:length(mod.fore)
         @test 0.85 < mod.back[i] / mod.fore[i] < 1.15
      end
   end

   # primer count model
   mod = JointBiasMod( 5 )
   test_uniform_bias!( mod )
   cnt = JointBiasCounter()
   @test cnt.count == 0.0
   hept = kmer_index_trailing(UInt16, dna"ATGAC")
   push!( cnt, hept )
   push!( cnt, hept )
   adjust!( cnt, mod )
   @test get(cnt) == value!( mod, hept )

   # gc content bias
   seq = dna"AAAAAGCGCG" ^ 5
   egc = ExpectedGC( seq )
   @test length(egc) == 20
   @test egc[ Int(div(gc_content(seq), 0.05)+1) ] == 1.0

end

@testset "Splice Graphs" begin
   gtf = IOBuffer("# gtf file test
chr0\tTEST\texon\t6\t20\t.\t+\t.\tgene_id \"one\"; transcript_id \"def\";
chr0\tTEST\texon\t31\t40\t.\t+\t.\tgene_id \"one\"; transcript_id \"def\";
chr0\tTEST\texon\t54\t62\t.\t+\t.\tgene_id \"one\"; transcript_id \"def\";
chr0\tTEST\texon\t76\t85\t.\t+\t.\tgene_id \"one\"; transcript_id \"def\";
chr0\tTEST\texon\t6\t40\t.\t+\t.\tgene_id \"one\"; transcript_id \"int1_alt3\";
chr0\tTEST\texon\t51\t62\t.\t+\t.\tgene_id \"one\"; transcript_id \"int1_alt3\";
chr0\tTEST\texon\t76\t85\t.\t+\t.\tgene_id \"one\"; transcript_id \"int1_alt3\";
chr0\tTEST\texon\t6\t20\t.\t+\t.\tgene_id \"one\"; transcript_id \"apa_alt5\";
chr0\tTEST\texon\t31\t40\t.\t+\t.\tgene_id \"one\"; transcript_id \"apa_alt5\";
chr0\tTEST\texon\t54\t65\t.\t+\t.\tgene_id \"one\"; transcript_id \"apa_alt5\";
chr0\tTEST\texon\t76\t90\t.\t+\t.\tgene_id \"one\"; transcript_id \"apa_alt5\";
chr0\tTEST\texon\t11\t20\t.\t+\t.\tgene_id \"single\"; transcript_id \"ex1_single\";
chr0\tTEST\texon\t11\t20\t.\t-\t.\tgene_id \"single_rev\"; transcript_id \"single_rev\";
chr0\tTEST\texon\t11\t20\t.\t-\t.\tgene_id \"kissing\"; transcript_id \"def_kiss\";
chr0\tTEST\texon\t21\t30\t.\t-\t.\tgene_id \"kissing\"; transcript_id \"def_kiss\";
chr0\tTEST\texon\t11\t30\t.\t-\t.\tgene_id \"kissing\"; transcript_id \"ret_kiss\";
")

#=   flat = IOBuffer("# refflat file test (gtfToGenePred -genePredExt test.gtf test.flat)
def\tchr0\t+\t5\t85\t85\t85\t4\t5,30,53,75,\t20,40,62,85,\t0\tone\tnone\tnone\t-1,-1,-1,-1,
int1_alt3\tchr0\t+\t5\t85\t85\t85\t3\t5,50,75,\t40,62,85,\t0\tone\tnone\tnone\t-1,-1,-1,
apa_alt5\tchr0\t+\t5\t90\t90\t90\t4\t5,30,53,75,\t20,40,65,90,\t0\tone\tnone\tnone\t-1,-1,-1,-1,
ex1_single\tchr0\t+\t10\t20\t10\t20\t1\t10,\t20,\t0\tsingle\tnone\tnone\t-1,
")=#

   gtfref  = load_gtf( gtf )
   #flatref = load_refflat( flat )

   @testset "Gene Annotations" begin

      @test gtfref["one"].don == (20,40,62,65)
      @test gtfref["one"].acc == (31,51,54,76)
      @test gtfref["one"].noveldon == tuple()
      @test gtfref["one"].novelacc == tuple()

   end

   gtf.ptr = 1
   bamreadr = open(BAM.Reader, "test.sort.bam", index="test.sort.bam.bai" )
   supgtfref = load_gtf( gtf, usebam=true, bamreader=Nullable(bamreadr), bamreads=2 )

   @testset "BAM Supplemented Annotation" begin

      @test supgtfref["one"].don == (20,40,45,62,65)
      @test supgtfref["one"].acc == (31,51,54,76)
      @test supgtfref["one"].noveldon == tuple(45)
      @test supgtfref["one"].novelacc == tuple()

   end

                               # fwd     rev
   buffer1   = dna"AAAAA"      # 1-5     96-100
   utr5      = dna"TTATT"      # 6-10    91-95
   exon1     = dna"GCGGATTACA" # 11-20   81-90
   int1      = dna"TTTTTTTTTT" # 21-30   71-80
   exon2     = dna"GCATTAGAAG" # 31-40   61-70
   int2      = dna"GGGGGGGGGG" # 41-50   51-60
   exon3alt3 = dna"CCT"        # 51-53   48-50
   exon3def  = dna"CTATGCTAG"  # 54-62   39-47
   exon3alt5 = dna"TTC"        # 63-65   36-38
   int3      = dna"CCCCCCCCCC" # 66-75   26-35
   exon4     = dna"TTAGACAAGA" # 76-85   16-25
   apa       = dna"AATAA"      # 86-90   11-15
   buffer2   = dna"AAAAAAAAAA" # 91-100  1-10

   fwd = buffer1 * utr5 * exon1 * int1 *
         exon2 * int2 *
         exon3alt3 * exon3def * exon3alt5 * int3 *
         exon4 * apa * buffer2
   rev = reverse_complement(fwd)

   genome = fwd * rev

   expected_one = utr5 * exon1 *
               int1 *
               exon2 *
               exon3alt3 *
               exon3def *
               exon3alt5 *
               exon4 *
               apa

   expected_sin = dna"GCGGATTACA"
   expected_kis = dna"AAAAAAAAAATGTAATCCGC"

   kmer_size = 2 # good test size

   graph_one = SpliceGraph( gtfref["one"], genome, kmer_size )
   graph_sin = SpliceGraph( gtfref["single"], genome, kmer_size )
   graph_rev = SpliceGraph( gtfref["single_rev"], genome, kmer_size )
   graph_kis = SpliceGraph( gtfref["kissing"], genome, kmer_size )

   @testset "Graph Building" begin
      @test graph_one.seq == expected_one
      @test graph_sin.seq == expected_sin
      @test graph_kis.seq == expected_kis

      @test length(graph_one.annopath) == length(gtfref["one"].reftx)
      @test graph_one.annopath[1] == IntSet([1,3,5,7]) # path of def
      @test graph_one.annopath[2] == IntSet([1,2,3,4,5,7]) # path of int1_alt3
      @test graph_one.annopath[3] == IntSet([1,3,5,6,7,8]) # path of apa_alt5
   end

   # Build Index (from index.jl)
   xcript  = dna""
   xoffset = Vector{UInt64}()
   xgenes  = Vector{GeneName}()
   xinfo   = Vector{GeneInfo}()
   xlength = Vector{Float64}()
   xgraph  = Vector{SpliceGraph}()

   runoffset = 0

   for g in keys(gtfref)
      println(STDERR, g)
      curgraph = SpliceGraph( gtfref[g], genome, kmer_size )
      xcript  *= curgraph.seq
      push!(xgraph, curgraph)
      push!(xgenes, g)
      push!(xinfo, gtfref[g].info )
      push!(xlength, gtfref[g].length )
      push!(xoffset, runoffset)
      runoffset += length(curgraph.seq)
   end

   fm = FMIndex(twobit_enc(xcript), 4, r=1, program=:SuffixArrays, mmap=true)

   edges = build_edges( xgraph, kmer_size )

   #println(edges.left)
   #println(edges.right)

   lib = GraphLib( xoffset, xgenes, xinfo, xlength, xgraph, edges, fm, true, kmer_size )

   @testset "Kmer Edges" begin
      left  = [dna"CA", dna"AG", dna"AG", dna"TC", dna"AA"]
      right = [dna"GC", dna"CC", dna"CT", dna"TT", dna"TG"]
      lkmer = map( x->kmer_index(SGKmer(x)), left )
      rkmer = map( x->kmer_index(SGKmer(x)), right )
      #println(lkmer[1])
      #println(rkmer[1])
      for i in 1:4^kmer_size
         if i in lkmer
            @test isassigned(edges.left, i)
            @test typeof(edges.left[i]) == Vector{SGNode}
            @test issorted(edges.left[i], lt=sortlt)
         else
            @test !isassigned(edges.left, i)
         end
         if i in rkmer
            @test isassigned(edges.right, i)
            @test typeof(edges.right[i]) == Vector{SGNode}
            @test issorted(edges.right[i], lt=sortlt)
         else
            @test !isassigned(edges.right, i)
         end
      end
      exon1_lind = lkmer[1]
      exon2_rind = rkmer[1]
      @test intersect( edges.left[exon1_lind], edges.right[exon2_rind] ) == edges.right[exon2_rind]
      @test edges.left[exon1_lind] ∩ edges.right[exon2_rind] == edges.right[exon2_rind]
   end

   @testset "Saving and Loading Index" begin
      println(STDERR, "Saving test index...")
      println(STDERR, lib)
      @timer open("test_index.jls", "w+") do io
         serialize( io, lib )
      end
      println(STDERR, "Loading test index...")
      @timer lib = open(deserialize, "test_index.jls")
      println(STDERR, lib)
   end

   # store general data structures outside of testset
   score_range = 0.05
   param = AlignParam( 1, 2, 4, 4, 4, 5, 1, 2, 1000, score_range, 0.7,
                          false, false, true, false, true )
   println(map( x->length(x.annoname), lib.graphs ))
   gquant = GraphLibQuant{SGAlignSingle,DefaultCounter}( lib )

   @testset "Alignment, SAM Format, Equivalence Classes" begin
      # reads
      fastq = IOBuffer("@1S10M%11,20%exon1
NGCGGATTACA
+
#BBBBBBBBBB
@1S9M%54,62%exon3def
NCTATGCTAG
+
#BBBBBBBBB
@1S15M%51,65%alt3-exon3-alt5
NCCTCTATGCTAGTTC
+
#BBBBBBBBBBBBBBB
@11M10N10M%10,40%exon1-exon2
NGCGGATTACAGCATTAGAAG
+
#BBBBBBBBBBBBBBBBBBBB
@5M10N6M%16,36%exon1trunc-exon2trunc
TTACAGCATTN
+
BBBBBBBBBB#
@10M33N9M1S%11,62%exon1-exon3def
GCGGATTACACTATGCTAGN
+
BBBBBBBBBBBBBBBBBBB#
@10M33N9M%11,62%exon1-exon3def:rc
CTAGCATAGTGTAATCCGC
+
BBBBBBBBBBBBBBBBBBB
@10M55N10M1S%11,85%exon1-exon4full
GCGGATTACATTAGACAAGAN
+
BBBBBBBBBBBBBBBBBBBB#
@11M55N4M%10,78%exon1-exon4_4bp
NGCGGATTACATTAG
+
#IIIIIIIIIIIIII
@2M55N10M%19,85%exon1_2bp-exon4:rc
TCTTGTCTAATG
+
IIIIIIIIIIII
")

      const DNASeqType = BioSequences.BioSequence{BioSequences.DNAAlphabet{2}}
      fqparse = FASTQ.Reader( BufferedInputStream(fastq), fill_ambiguous=DNA_C )
      reads  = allocate_fastq_records( 10 )
      read_chunk!( reads, fqparse )

      @test length(reads) == 10
      for r in reads
         fill!( r, 33 )
         println(STDERR, r)
         readname = string(r)

         align = ungapped_align( param, lib, r )
         println(STDERR, align)

         flush(STDERR)
         @test !isnull( align )
         @test length( align.value ) >= 1
         @test all(map( x->x.isvalid, align.value))
         scores = map( x->identity(x, length(r.sequence)), align.value )
         @test maximum(scores) - minimum(scores) <= score_range

         if length(search(readname, ":rc")) > 0
            @test align.value[1].strand == false
         else
            @test align.value[1].strand == true
         end

         best_ind = indmax(scores)

         ex_num = length(split(readname, '-', keep=false))
         @test length(align.value[best_ind].path) == ex_num

         count!( gquant, align.value[best_ind], 1.0 )

         curgene   = align.value[best_ind].path[1].gene
         firstnode = align.value[best_ind].path[1].node
         lastnode  = align.value[best_ind].path[end].node
         curgraph  = lib.graphs[ curgene ]

         cigar,positions = split(readname, '%', keep=false)[1:2]
         first,last   = split(positions, ',', keep=false) |> y->map(x->parse(Int,x), y)
         # Test SAM Format
         @test cigar == cigar_string( align.value[best_ind], curgraph, true, length(r.sequence) )[1]
         test_cigar,endpos = cigar_string( align.value[best_ind], curgraph, true, length(r.sequence) )
         #println(STDERR, "cigar = $test_cigar")
         @test first == Int(curgraph.nodecoord[firstnode] + (align.value[best_ind].offset - curgraph.nodeoffset[firstnode]))
         #@test last  == Int(curgraph.nodecoord[lastnode] + (endpos - curgraph.nodeoffset[lastnode]))
         # test readlength = number of M and S entries in cigar
         # test SAM offset is correct for both '+' and '-' genes.
         # test that cigar reversal works for '-' strand genes.
      end

      @testset "SGAlignContainer Hashing" begin

         @test SGAlignNode(1,1,zero(SGAlignScore)) == SGAlignNode(1,1,one(SGAlignScore))
         @test hash( SGAlignNode(1,1,zero(SGAlignScore)) ) == hash( SGAlignNode(1,1,one(SGAlignScore)) )

         @test SGAlignNode[SGAlignNode(1,1,zero(SGAlignScore))] == SGAlignNode[SGAlignNode(1,1,one(SGAlignScore))]
         @test hash( SGAlignNode[SGAlignNode(1,1,zero(SGAlignScore))] ) == hash( SGAlignNode[SGAlignNode(1,1,one(SGAlignScore))] )

         @test SGAlignSingle(SGAlignNode[SGAlignNode(1,1,zero(SGAlignScore))]) == SGAlignSingle(SGAlignNode[SGAlignNode(1,1,one(SGAlignScore))])
         @test hash( SGAlignSingle(SGAlignNode[SGAlignNode(1,1,zero(SGAlignScore))]) ) == hash( SGAlignSingle(SGAlignNode[SGAlignNode(1,1,one(SGAlignScore))]) )

         a = map( x->SGAlignSingle(SGAlignNode[SGAlignNode(x,x,zero(SGAlignScore))]), collect(1:5))
         ab = map( x->SGAlignSingle(SGAlignNode[SGAlignNode(x,x,one(SGAlignScore))]), collect(1:5))
         b = map( x->SGAlignPaired(SGAlignNode[SGAlignNode(x,x,zero(SGAlignScore))], SGAlignNode[SGAlignNode(x,x,zero(SGAlignScore))]), collect(1:5))
         bb = map( x->SGAlignPaired(SGAlignNode[SGAlignNode(x,x,one(SGAlignScore))], SGAlignNode[SGAlignNode(x,x,one(SGAlignScore))]), collect(1:5))

         da = Dict{SGAlignSingle,Int}()
         db = Dict{SGAlignPaired,Int}()

         function settwice!( d, a, ab )
            for i in a
               d[i] = 1
            end
            for i in ab
               d[i] = 2
            end
            for k in keys(d)
               @test d[k] == 2
            end
         end
         settwice!( da, a, ab )
         settwice!( db, b, bb )
         @test length(keys(da)) == length(a)
         @test length(keys(db)) == length(b)

      end

      @testset "SGAlignPath Overlap" begin

         single     = SGAlignNode[SGAlignNode(1,1,one(SGAlignScore))]
         single_two = SGAlignNode[SGAlignNode(1,2,one(SGAlignScore))]
         two_three  = SGAlignNode[SGAlignNode(1,2,one(SGAlignScore)), SGAlignNode(1,3,one(SGAlignScore))]
         larger     = SGAlignNode[SGAlignNode(1,1,one(SGAlignScore)), SGAlignNode(1,2,one(SGAlignScore)), SGAlignNode(1,3,one(SGAlignScore))]

         println(STDERR, @which single in single)

         @test single in single
         @test !(single in single_two)
         @test single_two in two_three
         @test !(single in two_three)
         @test single in larger
         @test single_two in larger
         @test two_three in larger

      end

      @testset "Isoform Equivalence Classes" begin

         is = IntSet([1, 2, 5, 6])
         @test in( 1, 2, is ) == true
         @test in( 1, 3, is ) == false
         @test in( 2, 3, is ) == false
         @test in( 2, 4, is ) == false
         @test in( 2, 5, is ) == true
         @test in( 5, 6, is ) == true
         @test in( 1, 6, is ) == false
         for i in [1,2,5,6]
            @test (i in is) == true
         end
         path = SGAlignNode[SGAlignNode(0x00000002, 0x00000001, SGAlignScore(0x02, 0x00, 0.0)), SGAlignNode(0x00000002, 0x00000007, SGAlignScore(0x0a, 0x00, 0.0))]
         @test (path in is) == false
         path = SGAlignNode[SGAlignNode(0x00000002, 0x00000001, SGAlignScore(0x02, 0x00, 0.0)), SGAlignNode(0x00000002, 0x00000002, SGAlignScore(0x0a, 0x00, 0.0))]
         @test (path in is) == true
         path = SGAlignNode[SGAlignNode(0x00000002, 0x00000001, SGAlignScore(0x02, 0x00, 0.0)), SGAlignNode(0x00000002, 0x00000002, SGAlignScore(0x0a, 0x00, 0.0)), SGAlignNode(0x00000002, 0x00000003, SGAlignScore(0x0a, 0x00, 0.0))]
         @test (path in is) == false
         path = SGAlignNode[SGAlignNode(0x00000002, 0x00000001, SGAlignScore(0x02, 0x00, 0.0)), SGAlignNode(0x00000002, 0x00000002, SGAlignScore(0x0a, 0x00, 0.0)), SGAlignNode(0x00000002, 0x00000005, SGAlignScore(0x0a, 0x00, 0.0))]
         @test (path in is) == true

         @test length(gquant.tpm) == 7
         @test length(gquant.count) == 7
         @test length(gquant.length) == 7
         @test gquant.geneidx == [0, 2, 5, 6]
         build_equivalence_classes!( gquant, lib, assign_long=true )
         println(STDERR, gquant)
         println(STDERR, gquant.classes)

      end

      @testset "MultiMapping Equivalence Classes" begin

         multi  = MultiMapping{SGAlignSingle,DefaultCounter}()
         aligns = SGAlignment[SGAlignment(0x0000000e, 0x01, SGAlignNode[SGAlignNode(0x00000002, 0x00000001, SGAlignScore(0x02, 0x00, 0.0)), SGAlignNode(0x00000002, 0x00000007, SGAlignScore(0x0a, 0x00, 0.0))], false, true),
                              SGAlignment(0x00000001, 0x01, SGAlignNode[SGAlignNode(0x00000004, 0x00000001, SGAlignScore(0x02, 0x00, 0.0))], false, true)]

         push!( multi, aligns, 1.0, gquant, lib )
         println(STDERR, multi.map)
         compats = collect(values(multi.map))
         @test compats[1].isdone == true
         @test sum(multi.iset) == 0

         multi  = MultiMapping{SGAlignSingle,DefaultCounter}()

         aligns = SGAlignment[SGAlignment(0x0000000e, 0x01, SGAlignNode[SGAlignNode(0x00000002, 0x00000001, SGAlignScore(0x02, 0x00, 0.0)), SGAlignNode(0x00000002, 0x00000002, SGAlignScore(0x0a, 0x00, 0.0))], false, true),
                              SGAlignment(0x00000001, 0x01, SGAlignNode[SGAlignNode(0x00000002, 0x00000001, SGAlignScore(0x02, 0x00, 0.0))], false, true)]

         push!( multi, aligns, 1.0, gquant, lib )
         println(STDERR, multi.map)
         compats = collect(values(multi.map))
         @test compats[1].isdone == false
         @test [3,4,5] == compats[1].class
         @test sum(multi.iset) == 0

         # Assignment
         compats[1].count = 10.0
         prev_node = get(gquant.quant[2].node[1])
         println(STDERR, "prev_quant = $(gquant.quant[2])")
         interv = Interval{ExonInt}( 1, 2 )
         prev_edge = get(get( gquant.quant[2].edge, interv, IntervalValue(0,0,zero(DefaultCounter)) ).value)

         println(gquant.quant[2].edge)

         assign_ambig!( gquant, lib, multi )
         println(gquant.quant[2].edge)

         println(STDERR, "cur_quant = $(gquant.quant[2])")
         @test get(gquant.quant[2].node[1]) == prev_node + 7.5
         @test get(get(gquant.quant[2].edge, interv, IntervalValue(0,0,default(DefaultCounter)) ).value) == prev_edge + 2.5
      end

      @testset "TPM" begin
         multi  = MultiMapping{SGAlignSingle,DefaultCounter}()
         calculate_tpm!( gquant, readlen=20 )
         set_gene_tpm!( gquant, lib )
         orig = deepcopy(gquant.tpm)
         iter = gene_em!( gquant, multi, sig=1, readlen=20, maxit=10000 )
         @test orig != gquant.tpm
         @test iter < 5000
         @test sum(gquant.tpm) == 1000000
      end

      @testset "Output Files" begin
         output_junctions( "test.jnc.gz", lib, gquant )
         output_stats( "test.map.gz", lib, gquant, param, "test_index.jls", 10, 10, 0, 20, "vTEST" )
         output_exons( "test_index.exons.gz", lib )
      end

      function parse_edge( str::S ) where S <: AbstractString
         s = split( str, ['-',':'] )
         (parse(Int, String(s[1])), parse(Int, String(s[2])), parse(Float64, String(s[3])))
      end

      @testset "Graph Enumeration" begin
         #= Real Test based on node 23 (here annotated as 1) in ENSMUSG00000038685
         inclusion edges:
         1-2, 2-3, 3-4, 4-5, 5-6, 6-7, 7-8, 8-9, 9-10
         exclusion edges:
         1-3, 2-4, 3-6, 3-10, 4-6, 8-10
         expected paths:
         =#
         expected_paths = [IntSet([1,2,3,4,5,6,7,8,9,10]) # (full)
         IntSet([1,3,4,5,6,7,8,9,10]) # (1-3)
         IntSet([1,2,4,5,6,7,8,9,10]) # (2-4)
         IntSet([1,2,3,6,7,8,9,10])   # (3-6)
         IntSet([1,3,6,7,8,9,10])     # (3-6)
         IntSet([1,2,3,10])           # (3-10)
         IntSet([1,3,10])             # (3-10)
         IntSet([1,2,3,4,6,7,8,9,10]) # (4-6)
         IntSet([1,3,4,6,7,8,9,10])   # (4-6)
         IntSet([1,2,4,6,7,8,9,10])   # (4-6)
         IntSet([1,2,3,4,5,6,7,8,10]) # (8-10)
         IntSet([1,3,4,5,6,7,8,10])   # (8-10)
         IntSet([1,2,4,5,6,7,8,10])   # (8-10)
         IntSet([1,2,3,6,7,8,10])     # (8-10)
         IntSet([1,3,6,7,8,10])       # (8-10)
         IntSet([1,2,3,4,6,7,8,10])   # (8-10)
         IntSet([1,3,4,6,7,8,10])     # (8-10)
         IntSet([1,2,4,6,7,8,10])]    # (8-10)
         edgestr = "1-2:616.0,1-3:4.0,2-3:569.0,2-4:20.0,3-4:629.0,3-6:1.0,3-10:3.0,4-5:1.0,4-6:664.0,5-6:5.0,6-7:789.0,7-8:790.0,8-9:2.0,8-10:606.0,9-10:15.0"
         edgespl  = split( edgestr, ',' )
         edges    = IntervalMap{Int,Float64}()
         psigraph = PsiGraph( Vector{Float64}(), Vector{Float64}(),
                              Vector{Float64}(), Vector{IntSet}(), 1, 10 )
         for s in split( edgestr, ',' )
            f,l,v = parse_edge( s )
            edges[(f,l)] = v
            push!( psigraph.psi, 0.0 )
            push!( psigraph.length, 1.0 )
            push!( psigraph.count, v )
            push!( psigraph.nodes, IntSet([f,l]) )
         end

         res = reduce_graph( psigraph )
         for p in expected_paths
            @test p in res.nodes
         end
         @test length(expected_paths) == length(res.nodes)
      end

      @testset "Event Building" begin
         out = IOBuffer(true,true)
         bs  = BufferedOutputStream(out)
         output_psi_header( bs )
         for g in 1:length(lib.graphs)
            name = lib.names[g]
            chr  = lib.info[g].name
            strand = lib.info[g].strand ? '+' : '-'
            _process_events( bs, lib.graphs[g], gquant.quant[g], (name,chr,strand), readlen=20, isnodeok=false )
         end
         flush(bs)
         seek(out,0)
         for l in eachline(out)
            spl = split( l, '\t' )
            @test length(spl) == 14
            println(STDERR,l)
         end
      end

      @testset "EBI Accessions & HTTP Streaming" begin
         ebi_res = ident_to_fastq_url("SRR1199010") # small single cell file
         @test ebi_res.success
         @test !ebi_res.paired

         println(STDERR, "Streaming fastq file from http://timbitz.github.io/data/SRR1199010.fastq.gz...")
         parser, response = make_http_fqparser( "http://timbitz.github.io/data/SRR1199010.fastq.gz" )
         println(STDERR, "HTTP Response: $response")
         reads  = allocate_fastq_records( 50 )
         cnt    = 0
         while length(reads) > 0
            read_http_chunk!( reads, parser, response )
            for r in reads
               fill!( r, 33 )
            end
            cnt += length(reads)
         end
         @test cnt == 128482 # correct number of reads in file
         #run(`julia ../bin/whippet-quant.jl --ebi SRR1199010`)
      end
   end
end


text:
using SafeTestsets: @safetestset

@time @safetestset "hankel" begin include("hankel.jl") end

ext: jl


text:
"""
Broadcast a message to all connected web pages to be interpetted by WebSocket listener. For example, in JavaScript:

var sock = new WebSocket('ws://'+window.location.host);
sock.onmessage = function( message ){
    var msg = JSON.parse(message.data);
    console.log(msg);
}
"""
function broadcast(msg::Dict)
    for (sid,s) in sessions
        c = s.client
        if ~c.is_closed
            write(c, json(msg))
        end
    end
end
broadcast(t,msg) = broadcast(Dict("type"=>t,"data"=>msg))

"""
Send a message to the specified connection to be interpetted by WebSocket listener. For example, in JavaScript:

var sock = new WebSocket('ws://'+window.location.host);
sock.onmessage = function( message ){
    var msg = JSON.parse(message.data);
    console.log(msg);
}
"""
function message(client::WebSocket,msg::Dict)
    if ~client.is_closed
        write(client, json(msg))
    end
end
message(client::WebSocket,t,msg) = message(client,Dict("type"=>t,"data"=>msg))
function message(id::Int,t,msg)
    for (sid,s) in sessions
        if isequal(id,s.client.id)
            message(s.client,Dict("type"=>t,"data"=>msg))
        end
    end
end

"""
Block Julia control flow until until callback["notify"](name) is called.
"""
function block(f::Function,name)
    conditions[name] = Condition()
    f()
    wait(conditions[name])
    delete!(conditions,name)
end

"""Add a JS library to the current page from a url."""
function add_library(url)
    name = basename(url)
    block(name) do
        broadcast("script","""
            var script = document.createElement("script");
            script.src = "$(url)";
            script.onload = Pages.notify("$(name)");
            document.head.appendChild(script);
        """)
    end
end

type Element
    id::String
    tag::String
    name::String
    attr::Dict{String,String}
    style::Dict{String,String}
    html::String
    text::String
    parent::String

    function Element(;id = "",tag = "div",name = "element",attr = Dict{String,String}(),style = Dict{String,String}(),html = "",text = "",parent = "body")
        new(id,tag,name,attr,style,html,text,parent)
    end
end

function assign(io::IO,element::Element)
    # ==========================================================================
    # Add attributes to element
    for prop in ["attr","style"]
        field = getfield(element,Symbol(prop))
        for key in keys(field)
            print(io,"""
                $(element.name).$(prop)("$(key)","$(field[key])");
            """)
        end
    end
    # ==========================================================================
    # Add html & text to element
    for prop in ["html","text"]
        field = getfield(element,Symbol(prop))
        if !isempty(field)
            print(io,"""
                $(element.name).$(prop)("$(field)");
            """)
        end
    end
    takebuf_string(io)
end

function add(io::IO,element::Element)
    # ==========================================================================
    # Get or create element
    print(io,"""
        var parent = d3.select("$(element.parent)")
        var $(element.name) = null;
    """)
    if !isempty(element.id)
        print(io,"""
            var check = document.getElementById("$(element.id)");
            if (check === null) {
                $(element.name) = parent.append("$(element.tag)").attr("id","$(element.id)");
            } else {
                $(element.name) = d3.select("#$(element.id)");
            };
        """)
    else
        print(io,"""
            $(element.name) = d3.select("$(element.parent)").append("$(element.tag)");
        """)
    end
    print(io,assign(io,element))
    takebuf_string(io)
end
function add(element::Element)
    Pages.broadcast("script",add(IOBuffer(),element))
end

function append(io::IO,element::Element;parent = """d3.select("body")""")
    if isempty(element.name)
        print(io,"""
            $(parent).append("$(element.tag)");
        """)
    else
        print(io,"""
            $(element.name) = $(parent).append("$(element.tag)");
        """)
    end
    print(io,assign(io,element))
    takebuf_string(io)
end

function remove(io::IO,tag;parent = """d3.select("body")""")
    print(io,"""
        $(parent).selectAll("$(tag)").remove();
    """)
    takebuf_string(io)
end
function remove(tag;parent = """d3.select("body")""")
    Pages.broadcast("script",remove(IOBuffer(),tag,parent=parent))
end

function add_select(io::IO,options,element::Element)
    element.tag == "select" || error("Element must have tag = select.")
    print(io,add(io,element))
    print(io,remove(io,"option",parent=element.name))
    # print(io,"""
    #     $(element.name).selectAll("option").remove();
    # """)
    for key in keys(options)
        print(io,"""
            $(element.name).append("option").attr("value","$(key)").text("$(options[key])");
        """)
    end
    takebuf_string(io)
end
function add_select(options,element::Element)
    Pages.broadcast("script",add_select(IOBuffer(),options,element))
end

function add_table(io::IO,df::DataFrame;table = Element(tag="table",name="table"),tr = Element(tag="tr",name="row"),th = Element(tag="th",name="header"),td = Element(tag="td",name="cell"))
    table.tag == "table" || error("Element must have tag = table.")
    tr.tag == "tr" || error("Element must have tag = tr.")
    th.tag == "th" || error("Element must have tag = th.")
    td.tag == "td" || error("Element must have tag = td.")
    print(io,add(io,table))
    print(io,remove(io,"tr",parent=table.name))
    # ==========================================================================
    # Add header
    print(io,"""
        var $(tr.name) = null;
    """)
    print(io,append(io,tr,parent=table.name))
    print(io,"""
        var $(th.name) = null;
    """)
    for name in names(df)
        th.html = string(name)
        print(io,append(io,th,parent=tr.name))
    end
    # ==========================================================================
    # Add data
    print(io,"""
        var $(td.name) = null;
    """)
    for irow in 1:size(df,1)
        row = df[irow,:]
        print(io,append(io,tr,parent=table.name))
        for name in names(df)
            td.html = string(row[name][1])
            print(io,append(io,td,parent=tr.name))
        end
    end
    takebuf_string(io)
end
function add_table(df::DataFrame;table = Element(tag="table",name="table"),tr = Element(tag="tr",name="row"),th = Element(tag="th",name="header"),td = Element(tag="td",name="cell"))
    Pages.broadcast("script",add_table(IOBuffer(),df,table=table,tr=tr,th=th,td=td))
end


text:
@testset "Multiple modifications" begin
    model = POI.Optimizer(MOI.Utilities.Model{Float64}())

    x = MOI.add_variables(model, 3)

    saf = MOI.ScalarAffineFunction(
        [
            MOI.ScalarAffineTerm(1.0, x[1]),
            MOI.ScalarAffineTerm(1.0, x[2]),
            MOI.ScalarAffineTerm(1.0, x[3]),
        ],
        0.0,
    )
    ci1 = MOI.add_constraint(model, saf, MOI.LessThan(1.0))
    ci2 = MOI.add_constraint(model, saf, MOI.LessThan(2.0))

    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        saf,
    )

    fc1 = MOI.get(model, MOI.ConstraintFunction(), ci1)
    @test MOI.coefficient.(fc1.terms) == [1.0, 1.0, 1.0]
    fc2 = MOI.get(model, MOI.ConstraintFunction(), ci2)
    @test MOI.coefficient.(fc2.terms) == [1.0, 1.0, 1.0]
    obj = MOI.get(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
    )
    @test MOI.coefficient.(obj.terms) == [1.0, 1.0, 1.0]

    changes_cis = [
        MOI.ScalarCoefficientChange(MOI.VariableIndex(1), 4.0)
        MOI.ScalarCoefficientChange(MOI.VariableIndex(1), 0.5)
        MOI.ScalarCoefficientChange(MOI.VariableIndex(3), 2.0)
    ]
    MOI.modify(model, [ci1, ci2, ci2], changes_cis)

    fc1 = MOI.get(model, MOI.ConstraintFunction(), ci1)
    @test MOI.coefficient.(fc1.terms) == [4.0, 1.0, 1.0]
    fc2 = MOI.get(model, MOI.ConstraintFunction(), ci2)
    @test MOI.coefficient.(fc2.terms) == [0.5, 1.0, 2.0]

    changes_obj = [
        MOI.ScalarCoefficientChange(MOI.VariableIndex(1), 4.0)
        MOI.ScalarCoefficientChange(MOI.VariableIndex(2), 10.0)
        MOI.ScalarCoefficientChange(MOI.VariableIndex(3), 2.0)
    ]
    MOI.modify(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        changes_obj,
    )

    obj = MOI.get(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
    )
    @test MOI.coefficient.(obj.terms) == [4.0, 10.0, 2.0]
end


text:
### A Pluto.jl notebook ###
# v0.16.0

using Markdown
using InteractiveUtils

# ╔═╡ 5072980b-06f4-472d-8801-786d9ea2848c
using PlutoUI, Plots, LinearAlgebra

# ╔═╡ 8b5a261a-9b2f-4d84-a3d8-ef34337de6d8
plotly()

# ╔═╡ f381f9e8-d165-4ae9-a0cb-2210e2faa453
TableOfContents(title="📚 Sadržaj", aside=true)

# ╔═╡ be21625e-1cf1-11eb-2efb-0b857e0bbd1e
md"
# Analiza trokuta

Zadan je trokut $\triangle ABC$. Kutove uz vrhove $A$, $B$ i $C$ označavamo s $\alpha$, $\beta$ i $\gamma$, redom.

Za zadani trokut automatski se računaju:

* opseg,
* površina,
* kutovi u stupnjevima,
* težišnice i težište,
* visine i sjecište visina,
* simetrale stranica, središte i radijus opisane kružnice, i
* simetrale kuteva, središte i radijus upisane kružnice,

te nacrtaju pripadne slike.

__Napomena__: Točke definiramo kao `Tuple()`, na primjer `T=(x,y,z)`. Radi računanja vektorskog produkta porebno je `Tuple()` konvertirati u tip `Vector{}` što se radi s naredbom `collect()`.
"

# ╔═╡ b40e846c-1cf2-11eb-0375-cd3f3e221358
# Zadajmo trokut ABC
begin
	A=(1,2,5)
	B=(-2,0,3)
	C=(2,1,-1)
end

# ╔═╡ 9bcc7de0-1d3e-11eb-1f26-49ae5603bd60
begin
	# Definicije za baratanje s točkama
	import Base.+
	+(A::Tuple, B::Tuple)=A.+B
	import Base.-
	-(A::Tuple, B::Tuple)=A.-B
	import Base./
	/(A::Tuple, B::Number)=A./B
	import LinearAlgebra.×
	×(A::Tuple,B::Tuple)=×(collect(A),collect(B))
	×(A::Array,B::Tuple)=×(A,collect(B))
end

# ╔═╡ 1bed9df2-1cf3-11eb-3da5-0ff167b35323
begin
	# Plots.markersizes=2
	mesh3d([A,B,C],legend=false, title="Trokut",color=:lightgray)
	scatter!([A,B,C],ms=2) # ms = markersize
	plot!(xlabel="x",ylabel="y",zlabel="z",aspect_ratio=:equal)
end

# ╔═╡ 52805212-1cf6-11eb-148c-b118898c3bd1
md"
## Opseg

__Opseg__ je zbroj duljina stranica.
"

# ╔═╡ 5e303fa0-1cf6-11eb-24f7-3b5da0b2957d
function opseg(A,B,C)
	return norm(B-A)+norm(C-A)+norm(C-B)
end

# ╔═╡ 7ce59f58-1cf6-11eb-138d-2fd2931a68ac
O=opseg(A,B,C)

# ╔═╡ 82a6f1a8-1cf6-11eb-0ead-71c9821841fe
md"
## Površina

$P=\frac{1}{2} | \overrightarrow{AB}\times \overrightarrow{AC}|$
"

# ╔═╡ 8c8bedb8-1cf6-11eb-28e0-d93fec161bb0
function površina(A,B,C)
	return norm((B-A)×(C-A))/2
end

# ╔═╡ b6a484e8-1cf6-11eb-2824-619bc952261c
P=površina(A,B,C)

# ╔═╡ 526b20c0-1d25-11eb-3167-934f431a9db8
md"
## Kutovi

$\alpha=\arccos \frac{\overrightarrow{AB}\cdot\overrightarrow{AC}}{|\overrightarrow{AB}||\overrightarrow{AC}|},\quad
\beta=\arccos \frac{\overrightarrow{BA}\cdot\overrightarrow{BC}}{|\overrightarrow{BA}||\overrightarrow{BC}|},\quad
\gamma=\arccos \frac{\overrightarrow{CA}\cdot\overrightarrow{CB}}{|\overrightarrow{CA}||\overrightarrow{CB}|}$
"

# ╔═╡ 229190e0-1d26-11eb-06ad-f31d46ff5432
function kutovi(A,B,C)
	# Funkcija računa kuteve u stupnjevima
	α=acosd((B-A)⋅(C-A)/(norm(B-A)*norm(C-A)))
	β=acosd((A-B)⋅(C-B)/(norm(A-B)*norm(C-B)))
	γ=acosd((A-C)⋅(B-C)/(norm(A-C)*norm(B-C)))
	return α,β,γ
end

# ╔═╡ 4719fba0-1d26-11eb-2267-b361bf7918e9
α,β,γ=kutovi(A,B,C)

# ╔═╡ 9d83b7b0-1d26-11eb-1210-f5da77878fde
# Provjera
α+β+γ

# ╔═╡ cabb2658-1cf6-11eb-2ffd-2dc8e7df8a46
md"
## Težište

__Težište__ je presjek težišnica, a __težišnica__ je spojnica vrha sa središtem nasuprotne stranice. Središte stranice nasuprot vrhu $A$ označavamo s $A_1$, itd.

Trebamo definirati strukturu `pravac` kao i funkciju koja računa presjek dvaju pravaca.
"

# ╔═╡ f8887980-1d43-11eb-248e-43e8c5a15c57
struct pravac
	T::Tuple # Točka
    s::Vector # vektor smjera
end

# ╔═╡ 549c4080-1d44-11eb-0ca9-cb31ea57b04c
function presjek(p::pravac,q::pravac)
	t=[p.s -q.s]\collect(q.T-p.T)
	return p.T+t[1].*Tuple(p.s)
end

# ╔═╡ 30cb1f98-1cf7-11eb-08d4-9b12985dacfd
function težište(A,B,C)
	A₁=(B+C)/2
	B₁=(A+C)/2
	C₁=(A+B)/2
	# Definirajmo težišnice iz vrhova A i B
	tA=pravac(A,collect(A₁-A))
	tB=pravac(B,collect(B₁-B))
	# Težište je presjek težišnica
	T=presjek(tA,tB)
	return A₁,B₁,C₁,T
end

# ╔═╡ 8590bc18-1cf7-11eb-302d-155630460d84
A₁,B₁,C₁,T=težište(A,B,C)

# ╔═╡ 996e3846-1cf7-11eb-11dc-dfd833157463
# Nacrtajmo polovišta nasuprotnih stranica i spojnice, te težište
begin
	mesh3d([A,B,C],legend=false, title="Težište",color=:lightgray)
	scatter!([A₁,B₁,C₁],ms=2)
	plot!([A,A₁])
	plot!([B,B₁])
	plot!([C,C₁],aspect_ratio=:equal)
	scatter!(T,color=:black,ms=2)
end

# ╔═╡ 6095a900-1d42-11eb-0f78-c79765cb47dc
md"
## Sjecište visina

__Sjecište visina__ ili __ortocentar__ je sjecište visina, a __visina__ je okomica iz vrha na suprotnu stranicu.
"

# ╔═╡ a56e2c00-1d42-11eb-292c-55752ce03ea6
function ortocentar(A,B,C)
	# Visina iz vrha A
	vA=pravac(A,((C-B)×(A-B))×(C-B))
	# Pravac kroz točke B i C
	pBC=pravac(B,collect(C-B))
	# Točka A₁ u kojoj visina iz vrha A siječe nasuprotnu stranicu
	A₁=presjek(vA,pBC)

	# Visina iz vrha B
	vB=pravac(B,((C-A)×(C-B))×(C-A))
	# Pravac kroz točke A i C
	pAC=pravac(A,collect(C-A))
	# Točka B₁ u kojoj visina iz vrha A siječe nasuprotnu stranicu
	B₁=presjek(vB,pAC)

	# Visina iz vrha C
	vC=pravac(C,((B-A)×(C-A))×(B-A))
	# Pravac kroz točke A i B
	pAB=pravac(A,collect(B-A))
	# Točka C₁ u kojoj visina iz vrha A siječe nasuprotnu stranicu
	C₁=presjek(vC,pAB)

	# Ortocentar
	O=presjek(vA,vB)
	return A₁,B₁,C₁,O
end

# ╔═╡ 85db7380-1d46-11eb-215a-41b96fe29cee
A₂,B₂,C₂,Oc=ortocentar(A,B,C)

# ╔═╡ 8e7d8590-1d42-11eb-031e-6d3fa68ae707
begin
	mesh3d([A,B,C],legend=false, title="Sjecište visina",color=:lightgray )
	scatter!([A,B,C],ms=2)
	scatter!([A₂,B₂,C₂],ms=2)
	scatter!(Oc,color=:black,ms=2)
	plot!(xlabel="x",ylabel="y",zlabel="z")
	# Provjera sijeku li se visine u istoj točki
	plot!([A,A₂])
	plot!([B,B₂])
	plot!([C,C₂])
	# Provjera ako je ortocentar izvan trokuta
	plot!([A,Oc])
	plot!([B,Oc])
	plot!([C,Oc],aspect_ratio=:equal)
end

# ╔═╡ dd7c1700-1d48-11eb-3cd8-07985f0d50c5
md"
## Opisana kružnica

Središte $S$ __opisane kružnice__ je sjecište __simetrala stranica__, a radijus $r$ je udaljenost od središta do bilo kojeg vrha.

Vektor smjera simetrale stranice se računa pomoću vektorsko-vektorskog produkta.
"

# ╔═╡ 184d7680-1d49-11eb-2505-b95fc1316ea0
function opisana_kružnica(A,B,C)
	# Simetrala stranice BC
	sBC=pravac((B+C)/2,((C-B)×(A-B))×(C-B))
	# Simetrala stranice AC
	sAC=pravac((A+C)/2,((C-A)×(C-B))×(C-A))
	# Simetrala stranice AB
	sAB=pravac((A+B)/2,((B-A)×(C-A))×(B-A))
	# Središte opisane kružnice
	S=presjek(sBC,sAC)
	# Radijus opisane kružnice
	r=norm(S-A)
	return S,r
end

# ╔═╡ dce70d80-1d49-11eb-2287-097c591556b3
S,r=opisana_kružnica(A,B,C)

# ╔═╡ 3066e240-1e09-11eb-3960-cf7fcafd4afb
md"
### Kružnica u prostoru

Točke na kružnici sa središtem u točki $S=(s_x,s_y,s_z)$ radijusa $r$ koja leži u ravnini $\mathcal{R}$ računaju su formulom

$$
\begin{pmatrix}x \\ y \\z \end{pmatrix}=\begin{pmatrix}s_x \\ s_y \\s_z \end{pmatrix} + (r\cos\phi)  \vec{u}_0 +(r \sin\phi )\vec{v}_0, \quad \phi\in[0,2\pi],$$

gdje su $\vec{u}_0$ i $\vec{v}_0$ međusobno okomiti jedinični vektori paralelni s ravninom $\mathcal{R}$.
"

# ╔═╡ 7cecba42-1e09-11eb-2c71-19ca6bb64aa1
function kružnica(A,B,C,S,r)
	u=collect(B-A)
	v=((B-A)×(C-A))×(B-A)
	u₀=u/norm(u)
	v₀=v/norm(v)
	Φ=range(0,stop=2*pi,length=201)
	return [S+r*cos(ϕ).*Tuple(u₀)+r*sin(ϕ).*Tuple(v₀) for ϕ in Φ]
end

# ╔═╡ ea65cb90-1d49-11eb-2743-0d5882b43328
begin
	mesh3d([A,B,C],legend=false, title="Opisana kružnica",color=:lightgray )
	scatter!([A,B,C],ms=2)
	scatter!(S,color=:black,ms=2)
	scatter!([(A+B)/2,(A+C)/2,(B+C)/2],ms=2)
	plot!(xlabel="x",ylabel="y",zlabel="z")
	# Vizualna provjera sijeku li se simetrale u istoj točki
	plot!([(A+B)/2,S])
	plot!([(A+C)/2,S])
	plot!([(B+C)/2,S])
	# Crtanje kružnice
    Cir=kružnica(A,B,C,S,r)
	plot!(Cir,aspect_ratio=:equal)
end

# ╔═╡ 8ff034e0-1df1-11eb-35b4-9b9e3178bef9
md"
## Upisana kružnica

Središte $S$ __upisane kružnice__ je sjecište __simetrala kutova__, a radijus $r$ je udaljenost od središta do bilo koje stranice.

Vektor smjera simetrale kuta se računa kao zbroj jediničnih vektora priležnih stranica (dijagonala romba raspolavlja kut).
"

# ╔═╡ af840200-1df1-11eb-2079-6f31dab59923
function upisana_kružnica(A,B,C)
	# Simetrala kuta α
	sα=pravac(A,collect((B-A)/norm(B-A)+(C-A)/norm(C-A)))
	# Simetrala kuta β
	sβ=pravac(B,collect((A-B)/norm(A-B)+(C-B)/norm(C-B)))
	# Simetrala kuta γ
	sγ=pravac(C,collect((C-A)/norm(C-A)+(C-B)/norm(C-B)))

	# Središte upisane kružnice
	S=presjek(sα,sβ)
	# Radijus upisane kružnice
	# Okomica na stranicu AB kroz središte S
	oAB=pravac(S,((B-A)×(C-A))×(B-A))
	# Pravac kroz stranicu AB
	pAB=pravac(A,collect(B-A))
	# Sjecište okomice i stranice AB
	T=presjek(oAB,pAB)
	r=norm(S-T)
	return S,r
end

# ╔═╡ 37382cf0-1dfa-11eb-0ff7-e16314375ea7
Su,ru=upisana_kružnica(A,B,C)

# ╔═╡ bc84dd80-1df1-11eb-30fd-45ff1ebe64d1
begin
	mesh3d([A,B,C],legend=false, title="Upisana kružnica",color=:lightgray )
	scatter!([A,B,C],ms=2)
	scatter!(Su,color=:black,ms=2)
	plot!(xlabel="x",ylabel="y",zlabel="z")
	# Vizualna provjera sijeku li se simetrale u istoj točki
	plot!([A,Su])
	plot!([B,Su])
	plot!([C,Su])
	# Crtanje kružnice
	Cir₁=kružnica(A,B,C,Su,ru)
	plot!(Cir₁,aspect_ratio=:equal)
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
Plots = "~1.22.3"
PlutoUI = "~0.7.14"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[Adapt]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "84918055d15b3114ede17ac6a7182f68870c16f7"
uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
version = "3.3.1"

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[Bzip2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "19a35467a82e236ff51bc17a3a44b69ef35185a2"
uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0"
version = "1.0.8+0"

[[Cairo_jll]]
deps = ["Artifacts", "Bzip2_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"]
git-tree-sha1 = "f2202b55d816427cd385a9a4f3ffb226bee80f99"
uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a"
version = "1.16.1+0"

[[ColorSchemes]]
deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random"]
git-tree-sha1 = "a851fec56cb73cfdf43762999ec72eff5b86882a"
uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
version = "3.15.0"

[[ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
version = "0.12.8"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "31d0151f5716b655421d9d75b7fa74cc4e744df2"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.39.0"

[[CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[Contour]]
deps = ["StaticArrays"]
git-tree-sha1 = "9f02045d934dc030edad45944ea80dbd1f0ebea7"
uuid = "d38c429a-6771-53c6-b99e-75d170b6e991"
version = "0.5.7"

[[DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "7d9d316f04214f7efdbb6398d545446e246eff02"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.10"

[[DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[EarCut_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "3f3a2501fa7236e9b911e0f7a588c657e822bb6d"
uuid = "5ae413db-bbd1-5e63-b57d-d24a61df00f5"
version = "2.2.3+0"

[[Expat_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b3bfd02e98aedfa5cf885665493c5598c350cd2f"
uuid = "2e619515-83b5-522b-bb60-26c02a35a201"
version = "2.2.10+0"

[[FFMPEG]]
deps = ["FFMPEG_jll"]
git-tree-sha1 = "b57e3acbe22f8484b4b5ff66a7499717fe1a9cc8"
uuid = "c87230d0-a227-11e9-1b43-d7ebe4e7570a"
version = "0.4.1"

[[FFMPEG_jll]]
deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "LAME_jll", "Libdl", "Ogg_jll", "OpenSSL_jll", "Opus_jll", "Pkg", "Zlib_jll", "libass_jll", "libfdk_aac_jll", "libvorbis_jll", "x264_jll", "x265_jll"]
git-tree-sha1 = "d8a578692e3077ac998b50c0217dfd67f21d1e5f"
uuid = "b22a6f82-2f65-5046-a5b2-351ab43fb4e5"
version = "4.4.0+0"

[[FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[Fontconfig_jll]]
deps = ["Artifacts", "Bzip2_jll", "Expat_jll", "FreeType2_jll", "JLLWrappers", "Libdl", "Libuuid_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "21efd19106a55620a188615da6d3d06cd7f6ee03"
uuid = "a3f928ae-7b40-5064-980b-68af3947d34b"
version = "2.13.93+0"

[[Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[FreeType2_jll]]
deps = ["Artifacts", "Bzip2_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "87eb71354d8ec1a96d4a7636bd57a7347dde3ef9"
uuid = "d7e528f0-a631-5988-bf34-fe36492bcfd7"
version = "2.10.4+0"

[[FriBidi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "aa31987c2ba8704e23c6c8ba8a4f769d5d7e4f91"
uuid = "559328eb-81f9-559d-9380-de523a88c83c"
version = "1.0.10+0"

[[GLFW_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libglvnd_jll", "Pkg", "Xorg_libXcursor_jll", "Xorg_libXi_jll", "Xorg_libXinerama_jll", "Xorg_libXrandr_jll"]
git-tree-sha1 = "dba1e8614e98949abfa60480b13653813d8f0157"
uuid = "0656b61e-2033-5cc2-a64a-77c0f6c09b89"
version = "3.3.5+0"

[[GR]]
deps = ["Base64", "DelimitedFiles", "GR_jll", "HTTP", "JSON", "Libdl", "LinearAlgebra", "Pkg", "Printf", "Random", "Serialization", "Sockets", "Test", "UUIDs"]
git-tree-sha1 = "c2178cfbc0a5a552e16d097fae508f2024de61a3"
uuid = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
version = "0.59.0"

[[GR_jll]]
deps = ["Artifacts", "Bzip2_jll", "Cairo_jll", "FFMPEG_jll", "Fontconfig_jll", "GLFW_jll", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Libtiff_jll", "Pixman_jll", "Pkg", "Qt5Base_jll", "Zlib_jll", "libpng_jll"]
git-tree-sha1 = "ef49a187604f865f4708c90e3f431890724e9012"
uuid = "d2c73de3-f751-5644-a686-071e5b155ba9"
version = "0.59.0+0"

[[GeometryBasics]]
deps = ["EarCut_jll", "IterTools", "LinearAlgebra", "StaticArrays", "StructArrays", "Tables"]
git-tree-sha1 = "58bcdf5ebc057b085e58d95c138725628dd7453c"
uuid = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
version = "0.4.1"

[[Gettext_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "XML2_jll"]
git-tree-sha1 = "9b02998aba7bf074d14de89f9d37ca24a1a0b046"
uuid = "78b55507-aeef-58d4-861c-77aaff3498b1"
version = "0.21.0+0"

[[Glib_jll]]
deps = ["Artifacts", "Gettext_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "7bf67e9a481712b3dbe9cb3dac852dc4b1162e02"
uuid = "7746bdde-850d-59dc-9ae8-88ece973131d"
version = "2.68.3+0"

[[Graphite2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "344bf40dcab1073aca04aa0df4fb092f920e4011"
uuid = "3b182d85-2403-5c21-9c21-1e1f0cc25472"
version = "1.3.14+0"

[[Grisu]]
git-tree-sha1 = "53bb909d1151e57e2484c3d1b53e19552b887fb2"
uuid = "42e2da0e-8278-4e71-bc24-59509adca0fe"
version = "1.0.2"

[[HTTP]]
deps = ["Base64", "Dates", "IniFile", "Logging", "MbedTLS", "NetworkOptions", "Sockets", "URIs"]
git-tree-sha1 = "14eece7a3308b4d8be910e265c724a6ba51a9798"
uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
version = "0.9.16"

[[HarfBuzz_jll]]
deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "Graphite2_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg"]
git-tree-sha1 = "8a954fed8ac097d5be04921d595f741115c1b2ad"
uuid = "2e76f6c2-a576-52d4-95c1-20adfe4de566"
version = "2.8.1+0"

[[HypertextLiteral]]
git-tree-sha1 = "72053798e1be56026b81d4e2682dbe58922e5ec9"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.0"

[[IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[IniFile]]
deps = ["Test"]
git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8"
uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
version = "0.5.0"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[IterTools]]
git-tree-sha1 = "05110a2ab1fc5f932622ffea2a003221f4782c18"
uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
version = "1.3.0"

[[IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[JpegTurbo_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "d735490ac75c5cb9f1b00d8b5509c11984dc6943"
uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8"
version = "2.1.0+0"

[[LAME_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "f6250b16881adf048549549fba48b1161acdac8c"
uuid = "c1c5ebd0-6772-5130-a774-d5fcae4a789d"
version = "3.100.1+0"

[[LZO_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "e5b909bcf985c5e2605737d2ce278ed791b89be6"
uuid = "dd4b983a-f0e5-5f8d-a1b7-129d4a5fb1ac"
version = "2.10.1+0"

[[LaTeXStrings]]
git-tree-sha1 = "c7f1c695e06c01b95a67f0cd1d34994f3e7db104"
uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
version = "1.2.1"

[[Latexify]]
deps = ["Formatting", "InteractiveUtils", "LaTeXStrings", "MacroTools", "Markdown", "Printf", "Requires"]
git-tree-sha1 = "a4b12a1bd2ebade87891ab7e36fdbce582301a92"
uuid = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
version = "0.15.6"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[Libffi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "761a393aeccd6aa92ec3515e428c26bf99575b3b"
uuid = "e9f186c6-92d2-5b65-8a66-fee21dc1b490"
version = "3.2.2+0"

[[Libgcrypt_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgpg_error_jll", "Pkg"]
git-tree-sha1 = "64613c82a59c120435c067c2b809fc61cf5166ae"
uuid = "d4300ac3-e22c-5743-9152-c294e39db1e4"
version = "1.8.7+0"

[[Libglvnd_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll", "Xorg_libXext_jll"]
git-tree-sha1 = "7739f837d6447403596a75d19ed01fd08d6f56bf"
uuid = "7e76a0d4-f3c7-5321-8279-8d96eeed0f29"
version = "1.3.0+3"

[[Libgpg_error_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "c333716e46366857753e273ce6a69ee0945a6db9"
uuid = "7add5ba3-2f88-524e-9cd5-f83b8a55f7b8"
version = "1.42.0+0"

[[Libiconv_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "42b62845d70a619f063a7da093d995ec8e15e778"
uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
version = "1.16.1+1"

[[Libmount_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "9c30530bf0effd46e15e0fdcf2b8636e78cbbd73"
uuid = "4b2f31a3-9ecc-558c-b454-b3730dcb73e9"
version = "2.35.0+0"

[[Libtiff_jll]]
deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Pkg", "Zlib_jll", "Zstd_jll"]
git-tree-sha1 = "340e257aada13f95f98ee352d316c3bed37c8ab9"
uuid = "89763e89-9b03-5906-acba-b20f662cd828"
version = "4.3.0+0"

[[Libuuid_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "7f3efec06033682db852f8b3bc3c1d2b0a0ab066"
uuid = "38a345b3-de98-5d2b-a5d3-14cd9215e700"
version = "2.36.0+0"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "5a5bc6bf062f0f95e62d0fe0a2d99699fed82dd9"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.8"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS]]
deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"]
git-tree-sha1 = "1c38e51c3d08ef2278062ebceade0e46cefc96fe"
uuid = "739be429-bea8-5141-9913-cc70e7f3736d"
version = "1.0.3"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Measures]]
git-tree-sha1 = "e498ddeee6f9fdb4551ce855a46f54dbd900245f"
uuid = "442fdcdd-2543-5da2-b0f3-8c86c306513e"
version = "0.3.1"

[[Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[NaNMath]]
git-tree-sha1 = "bfe47e760d60b82b66b61d2d44128b62e3a369fb"
uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
version = "0.3.5"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[Ogg_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "7937eda4681660b4d6aeeecc2f7e1c81c8ee4e2f"
uuid = "e7412a2a-1a6e-54c0-be00-318e2571c051"
version = "1.3.5+0"

[[OpenSSL_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "15003dcb7d8db3c6c857fda14891a539a8f2705a"
uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95"
version = "1.1.10+0"

[[Opus_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "51a08fb14ec28da2ec7a927c4337e4332c2a4720"
uuid = "91d4177d-7536-5919-b921-800302f37372"
version = "1.3.2+0"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[PCRE_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b2a7af664e098055a7529ad1a900ded962bca488"
uuid = "2f80f16e-611a-54ab-bc61-aa92de5b98fc"
version = "8.44.0+0"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "a8709b968a1ea6abc2dc1967cb1db6ac9a00dfb6"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.0.5"

[[Pixman_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b4f5d02549a10e20780a24fce72bea96b6329e29"
uuid = "30392449-352a-5448-841d-b1acce4e97dc"
version = "0.40.1+0"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[PlotThemes]]
deps = ["PlotUtils", "Requires", "Statistics"]
git-tree-sha1 = "a3a964ce9dc7898193536002a6dd892b1b5a6f1d"
uuid = "ccf2f8ad-2431-5c83-bf29-c5338b663b6a"
version = "2.0.1"

[[PlotUtils]]
deps = ["ColorSchemes", "Colors", "Dates", "Printf", "Random", "Reexport", "Statistics"]
git-tree-sha1 = "b084324b4af5a438cd63619fd006614b3b20b87b"
uuid = "995b91a9-d308-5afd-9ec6-746e21dbc043"
version = "1.0.15"

[[Plots]]
deps = ["Base64", "Contour", "Dates", "Downloads", "FFMPEG", "FixedPointNumbers", "GR", "GeometryBasics", "JSON", "Latexify", "LinearAlgebra", "Measures", "NaNMath", "PlotThemes", "PlotUtils", "Printf", "REPL", "Random", "RecipesBase", "RecipesPipeline", "Reexport", "Requires", "Scratch", "Showoff", "SparseArrays", "Statistics", "StatsBase", "UUIDs"]
git-tree-sha1 = "cfbd033def161db9494f86c5d18fbf874e09e514"
uuid = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
version = "1.22.3"

[[PlutoUI]]
deps = ["Base64", "Dates", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "d1fb76655a95bf6ea4348d7197b22e889a4375f4"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.14"

[[Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[Qt5Base_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Fontconfig_jll", "Glib_jll", "JLLWrappers", "Libdl", "Libglvnd_jll", "OpenSSL_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libxcb_jll", "Xorg_xcb_util_image_jll", "Xorg_xcb_util_keysyms_jll", "Xorg_xcb_util_renderutil_jll", "Xorg_xcb_util_wm_jll", "Zlib_jll", "xkbcommon_jll"]
git-tree-sha1 = "ad368663a5e20dbb8d6dc2fddeefe4dae0781ae8"
uuid = "ea2cea3b-5b76-57ae-a6ef-0a8af62496e1"
version = "5.15.3+0"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[RecipesBase]]
git-tree-sha1 = "44a75aa7a527910ee3d1751d1f0e4148698add9e"
uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
version = "1.1.2"

[[RecipesPipeline]]
deps = ["Dates", "NaNMath", "PlotUtils", "RecipesBase"]
git-tree-sha1 = "7ad0dfa8d03b7bcf8c597f59f5292801730c55b8"
uuid = "01d81517-befc-4cb6-b9ec-a95719d0359c"
version = "0.4.1"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.1.3"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Scratch]]
deps = ["Dates"]
git-tree-sha1 = "0b4b7f1393cff97c33891da2a0bf69c6ed241fda"
uuid = "6c6a2e73-6563-6170-7368-637461726353"
version = "1.1.0"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[Showoff]]
deps = ["Dates", "Grisu"]
git-tree-sha1 = "91eddf657aca81df9ae6ceb20b959ae5653ad1de"
uuid = "992d4aef-0814-514b-bc4d-f2e9a6c4116f"
version = "1.0.3"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "3c76dde64d03699e074ac02eb2e8ba8254d428da"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.2.13"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[StatsAPI]]
git-tree-sha1 = "1958272568dc176a1d881acb797beb909c785510"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.0.0"

[[StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "8cbbc098554648c84f79a463c9ff0fd277144b6c"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.33.10"

[[StructArrays]]
deps = ["Adapt", "DataAPI", "StaticArrays", "Tables"]
git-tree-sha1 = "2ce41e0d042c60ecd131e9fb7154a3bfadbf50d3"
uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
version = "0.6.3"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "fed34d0e71b91734bf0a7e10eb1bb05296ddbcd0"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.6.0"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[URIs]]
git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355"
uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
version = "1.3.0"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[Wayland_jll]]
deps = ["Artifacts", "Expat_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg", "XML2_jll"]
git-tree-sha1 = "3e61f0b86f90dacb0bc0e73a0c5a83f6a8636e23"
uuid = "a2964d1f-97da-50d4-b82a-358c7fce9d89"
version = "1.19.0+0"

[[Wayland_protocols_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Wayland_jll"]
git-tree-sha1 = "2839f1c1296940218e35df0bbb220f2a79686670"
uuid = "2381bf8a-dfd0-557d-9999-79630e7b1b91"
version = "1.18.0+4"

[[XML2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "1acf5bdf07aa0907e0a37d3718bb88d4b687b74a"
uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
version = "2.9.12+0"

[[XSLT_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgcrypt_jll", "Libgpg_error_jll", "Libiconv_jll", "Pkg", "XML2_jll", "Zlib_jll"]
git-tree-sha1 = "91844873c4085240b95e795f692c4cec4d805f8a"
uuid = "aed1982a-8fda-507f-9586-7b0439959a61"
version = "1.1.34+0"

[[Xorg_libX11_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxcb_jll", "Xorg_xtrans_jll"]
git-tree-sha1 = "5be649d550f3f4b95308bf0183b82e2582876527"
uuid = "4f6342f7-b3d2-589e-9d20-edeb45f2b2bc"
version = "1.6.9+4"

[[Xorg_libXau_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4e490d5c960c314f33885790ed410ff3a94ce67e"
uuid = "0c0b7dd1-d40b-584c-a123-a41640f87eec"
version = "1.0.9+4"

[[Xorg_libXcursor_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXfixes_jll", "Xorg_libXrender_jll"]
git-tree-sha1 = "12e0eb3bc634fa2080c1c37fccf56f7c22989afd"
uuid = "935fb764-8cf2-53bf-bb30-45bb1f8bf724"
version = "1.2.0+4"

[[Xorg_libXdmcp_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4fe47bd2247248125c428978740e18a681372dd4"
uuid = "a3789734-cfe1-5b06-b2d0-1dd0d9d62d05"
version = "1.1.3+4"

[[Xorg_libXext_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "b7c0aa8c376b31e4852b360222848637f481f8c3"
uuid = "1082639a-0dae-5f34-9b06-72781eeb8cb3"
version = "1.3.4+4"

[[Xorg_libXfixes_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "0e0dc7431e7a0587559f9294aeec269471c991a4"
uuid = "d091e8ba-531a-589c-9de9-94069b037ed8"
version = "5.0.3+4"

[[Xorg_libXi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll", "Xorg_libXfixes_jll"]
git-tree-sha1 = "89b52bc2160aadc84d707093930ef0bffa641246"
uuid = "a51aa0fd-4e3c-5386-b890-e753decda492"
version = "1.7.10+4"

[[Xorg_libXinerama_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll"]
git-tree-sha1 = "26be8b1c342929259317d8b9f7b53bf2bb73b123"
uuid = "d1454406-59df-5ea1-beac-c340f2130bc3"
version = "1.1.4+4"

[[Xorg_libXrandr_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll"]
git-tree-sha1 = "34cea83cb726fb58f325887bf0612c6b3fb17631"
uuid = "ec84b674-ba8e-5d96-8ba1-2a689ba10484"
version = "1.5.2+4"

[[Xorg_libXrender_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "19560f30fd49f4d4efbe7002a1037f8c43d43b96"
uuid = "ea2f1a96-1ddc-540d-b46f-429655e07cfa"
version = "0.9.10+4"

[[Xorg_libpthread_stubs_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "6783737e45d3c59a4a4c4091f5f88cdcf0908cbb"
uuid = "14d82f49-176c-5ed1-bb49-ad3f5cbd8c74"
version = "0.1.0+3"

[[Xorg_libxcb_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "XSLT_jll", "Xorg_libXau_jll", "Xorg_libXdmcp_jll", "Xorg_libpthread_stubs_jll"]
git-tree-sha1 = "daf17f441228e7a3833846cd048892861cff16d6"
uuid = "c7cfdc94-dc32-55de-ac96-5a1b8d977c5b"
version = "1.13.0+3"

[[Xorg_libxkbfile_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "926af861744212db0eb001d9e40b5d16292080b2"
uuid = "cc61e674-0454-545c-8b26-ed2c68acab7a"
version = "1.1.0+4"

[[Xorg_xcb_util_image_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "0fab0a40349ba1cba2c1da699243396ff8e94b97"
uuid = "12413925-8142-5f55-bb0e-6d7ca50bb09b"
version = "0.4.0+1"

[[Xorg_xcb_util_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxcb_jll"]
git-tree-sha1 = "e7fd7b2881fa2eaa72717420894d3938177862d1"
uuid = "2def613f-5ad1-5310-b15b-b15d46f528f5"
version = "0.4.0+1"

[[Xorg_xcb_util_keysyms_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "d1151e2c45a544f32441a567d1690e701ec89b00"
uuid = "975044d2-76e6-5fbe-bf08-97ce7c6574c7"
version = "0.4.0+1"

[[Xorg_xcb_util_renderutil_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "dfd7a8f38d4613b6a575253b3174dd991ca6183e"
uuid = "0d47668e-0667-5a69-a72c-f761630bfb7e"
version = "0.3.9+1"

[[Xorg_xcb_util_wm_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "e78d10aab01a4a154142c5006ed44fd9e8e31b67"
uuid = "c22f9ab0-d5fe-5066-847c-f4bb1cd4e361"
version = "0.4.1+1"

[[Xorg_xkbcomp_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxkbfile_jll"]
git-tree-sha1 = "4bcbf660f6c2e714f87e960a171b119d06ee163b"
uuid = "35661453-b289-5fab-8a00-3d9160c6a3a4"
version = "1.4.2+4"

[[Xorg_xkeyboard_config_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xkbcomp_jll"]
git-tree-sha1 = "5c8424f8a67c3f2209646d4425f3d415fee5931d"
uuid = "33bec58e-1273-512f-9401-5d533626f822"
version = "2.27.0+4"

[[Xorg_xtrans_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "79c31e7844f6ecf779705fbc12146eb190b7d845"
uuid = "c5fb5394-a638-5e4d-96e5-b29de1b5cf10"
version = "1.4.0+3"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[Zstd_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "cc4bf3fdde8b7e3e9fa0351bdeedba1cf3b7f6e6"
uuid = "3161d3a3-bdf6-5164-811a-617609db77b4"
version = "1.5.0+0"

[[libass_jll]]
deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "HarfBuzz_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "5982a94fcba20f02f42ace44b9894ee2b140fe47"
uuid = "0ac62f75-1d6f-5e53-bd7c-93b484bb37c0"
version = "0.15.1+0"

[[libfdk_aac_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "daacc84a041563f965be61859a36e17c4e4fcd55"
uuid = "f638f0a6-7fb0-5443-88ba-1cc74229b280"
version = "2.0.2+0"

[[libpng_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "94d180a6d2b5e55e447e2d27a29ed04fe79eb30c"
uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f"
version = "1.6.38+0"

[[libvorbis_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Ogg_jll", "Pkg"]
git-tree-sha1 = "c45f4e40e7aafe9d086379e5578947ec8b95a8fb"
uuid = "f27f6e37-5d2b-51aa-960f-b287f2bc3b7a"
version = "1.3.7+0"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"

[[x264_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4fea590b89e6ec504593146bf8b988b2c00922b2"
uuid = "1270edf5-f2f9-52d2-97e9-ab00b5d0237a"
version = "2021.5.5+0"

[[x265_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "ee567a171cce03570d77ad3a43e90218e38937a9"
uuid = "dfaa095f-4041-5dcd-9319-2fabd8486b76"
version = "3.5.0+0"

[[xkbcommon_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Wayland_jll", "Wayland_protocols_jll", "Xorg_libxcb_jll", "Xorg_xkeyboard_config_jll"]
git-tree-sha1 = "ece2350174195bb31de1a63bea3a41ae1aa593b6"
uuid = "d8fb68d0-12a3-5cfd-a85a-d49703b185fd"
version = "0.9.1+5"
"""

# ╔═╡ Cell order:
# ╠═5072980b-06f4-472d-8801-786d9ea2848c
# ╠═8b5a261a-9b2f-4d84-a3d8-ef34337de6d8
# ╠═f381f9e8-d165-4ae9-a0cb-2210e2faa453
# ╟─be21625e-1cf1-11eb-2efb-0b857e0bbd1e
# ╠═b40e846c-1cf2-11eb-0375-cd3f3e221358
# ╠═9bcc7de0-1d3e-11eb-1f26-49ae5603bd60
# ╠═1bed9df2-1cf3-11eb-3da5-0ff167b35323
# ╠═52805212-1cf6-11eb-148c-b118898c3bd1
# ╠═5e303fa0-1cf6-11eb-24f7-3b5da0b2957d
# ╠═7ce59f58-1cf6-11eb-138d-2fd2931a68ac
# ╟─82a6f1a8-1cf6-11eb-0ead-71c9821841fe
# ╠═8c8bedb8-1cf6-11eb-28e0-d93fec161bb0
# ╠═b6a484e8-1cf6-11eb-2824-619bc952261c
# ╟─526b20c0-1d25-11eb-3167-934f431a9db8
# ╠═229190e0-1d26-11eb-06ad-f31d46ff5432
# ╠═4719fba0-1d26-11eb-2267-b361bf7918e9
# ╠═9d83b7b0-1d26-11eb-1210-f5da77878fde
# ╟─cabb2658-1cf6-11eb-2ffd-2dc8e7df8a46
# ╠═f8887980-1d43-11eb-248e-43e8c5a15c57
# ╠═549c4080-1d44-11eb-0ca9-cb31ea57b04c
# ╠═30cb1f98-1cf7-11eb-08d4-9b12985dacfd
# ╠═8590bc18-1cf7-11eb-302d-155630460d84
# ╠═996e3846-1cf7-11eb-11dc-dfd833157463
# ╟─6095a900-1d42-11eb-0f78-c79765cb47dc
# ╠═a56e2c00-1d42-11eb-292c-55752ce03ea6
# ╠═85db7380-1d46-11eb-215a-41b96fe29cee
# ╠═8e7d8590-1d42-11eb-031e-6d3fa68ae707
# ╟─dd7c1700-1d48-11eb-3cd8-07985f0d50c5
# ╠═184d7680-1d49-11eb-2505-b95fc1316ea0
# ╠═dce70d80-1d49-11eb-2287-097c591556b3
# ╟─3066e240-1e09-11eb-3960-cf7fcafd4afb
# ╠═7cecba42-1e09-11eb-2c71-19ca6bb64aa1
# ╠═ea65cb90-1d49-11eb-2743-0d5882b43328
# ╟─8ff034e0-1df1-11eb-35b4-9b9e3178bef9
# ╠═af840200-1df1-11eb-2079-6f31dab59923
# ╠═37382cf0-1dfa-11eb-0ff7-e16314375ea7
# ╠═bc84dd80-1df1-11eb-30fd-45ff1ebe64d1
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_name: ivanslapnicar/Matematika

max_issues_repo_name: ivanslapnicar/Matematika

max_forks_repo_name: ivanslapnicar/Matematika


text:
using Dates

# The time unit is days (as floating point)
# Day 0 is taken at 1 February 2020
const BASE_DATE = Date(2020, 2, 1)
const BASE_DAYS = 0

#%% R₀ calculations
# Default values for R_0
const BaseR₀ = 2.3954

# Seasonal forcing parameter ϵ
const ϵ = Dict(:north => 0.2, :tropical => 0.0, :south => 0.2)

#-------------------------------------------------------------------------------------------------
#--
#-- EPIDEMIOLOGY
#-

#-- Susceptibility to contagion and transition from a compartment to another
const AgeGroup = ["0-9", "10-19", "20-29", "30-39", "40-49", "50-59", "60-69", "70-79", "80+"]
const zₐ =       [0.05,   0.05,   0.10,    0.15,    0.20,    0.25,    0.30,    0.40,    0.50]
const mₐ =       [0.01,   0.03,   0.03,    0.03,    0.06,    0.10,    0.25,    0.35,    0.50]
const cₐ =       [0.05,   0.10,   0.10,    0.15,    0.20,    0.25,    0.35,    0.45,    0.55]
const fₐ =       [0.30,   0.30,   0.30,    0.30,    0.30,    0.40,    0.40,    0.50,    0.50]

const nAgeGroup = length(AgeGroup)

#-- Transition times (all in days)
# Time to infectiousness (written t\_l)
const tₗ = 5.653

# Time to infectiousness (written t\_i)
const tᵢ = 2.676

# Time in hospital bed (not ICU)
const tₕ = 2.007

# Time in ICU
const tᵤ = 13.884

# Time in asymptomatic recovery
const tᵣ = 14

#-- Contagion multiplier
# asymptomatic individuals (compartment E)
const γₑ = 0.33407

# Infected / symptomatic individuals
const γᵢ=0.25587

# Severe symptoms
const γⱼ=0.47757

# Critical symptoms
const γₖ = 5.94806

# Asymptomatic recovery
const γᵣ = 0.2

#-- Fatality Multiplier
# In ICU
const δᵤ = 1.0

# In hospital
const δₗ = 0.82117

# Out of hospital
const δₖ = 5.54388

#-------------------------------------------------------------------------------------------------
# Create parameter vectors containing decribing each parameter
#            Name             Base value      Allowed range      Disease     +/- for optim
#                                                                specific?

const DEATH_AT_MODEL_START = 1.0

const DISEASE_BASE =
         [  ["r₀",            BaseR₀,         (2.0, 3.0),        1,          1.0],
            ["tₗ",            tₗ,             (1.0, 8.0),        1,          1.0],
            ["tᵢ",            tᵢ,             (1.0, 8.0),        1,          1.0],
            ["tₕ",            tₕ,             (1.0, 8.0),        1,          1.0],
            ["tᵤ",            tᵤ,             (1.0, 25.0),       1,          2.0],
            ["tᵣ",            tᵣ,             (0.01, 40.0),      1,          2.0],
            ["γₑ",            γₑ,             (0.0, 2.0),        1,          0.3],
            ["γᵢ",            γᵢ,             (0.0, 2.0),        1,          0.0],
            ["γⱼ",            γⱼ,             (1.0, 1.0),        1,          1.0],
            ["γₖ",            γₖ,             (0.0, 10.0),       1,          2.0],
            ["γᵣ",            γᵣ,             (0.0, 10.0),       1,          2.0],
            ["δₖ",            δₖ,             (0.5, 3.0),        1,          0.5],
            ["δₗ",            δₗ,             (0.5, 9.0),        1,          1.0],
            ["δᵤ",            δᵤ,             (1.0, 1.0),        1,          0.0]]

const DISEASE_N     = size(DISEASE_BASE)[1]
const DISEASE_NAMES = [DISEASE_BASE[i][1]                        for i in 1:DISEASE_N]
const DISEASE_INIT  = [DISEASE_BASE[i][2]                        for i in 1:DISEASE_N]
const DISEASE_RANGE = Tuple{Float64, Float64}[DISEASE_BASE[i][3] for i in 1:DISEASE_N]
const DISEASE_OPTIM = [DISEASE_BASE[i][5] for i in 1:DISEASE_N]

const COUNTRY_BASE =
            # modelStart shifts between model dates which always start at 0 and real
            # dates where the model starts at a different date
            # add modelStart to go from model time to actual (eg -30 means the model starts
            # 30 days before BASE_DATE)
            # substract to go from real (eg. date of deaths) to model time
            # The range of the value is approximated by approximateModelStartRange()
         [  ["modelStart",    0.0,            (-90.0, 45.0),     0,          5.0],
            ["infectedM",     25.0,           (1.0, 5000.0),     0,          10.0],
            ["infectiousM",   25.0,           (1.0, 5000.0),     0,          10.0],

            # Mitigation profile with 4 points
            ["mv0",           1.0,            (1.0, 1.0),        0,          0.0],
            ["mv1",           0.8,            (0.1, 1.5),        0,          0.3],
            ["mv2",           0.5,            (0.1, 1.5),        0,          0.3],
            ["mv3",           0.5,            (0.1, 1.5),        0,          0.3],
            ["mv4",           0.5,            (0.1, 1.5),        0,          0.3],
            ["mv5",           0.5,            (0.1, 1.5),        0,          0.3],
            ["mv6",           0.5,            (0.1, 1.5),        0,          0.3],
            ["mv7",           0.5,            (0.1, 1.5),        0,          0.3],
            ["mv8",           0.5,            (0.1, 1.5),        0,          0.3],
            ["mv9",           0.5,            (0.1, 1.5),        0,          0.3]]

const COUNTRY_N = size(COUNTRY_BASE)[1]
const COUNTRY_NAMES = [COUNTRY_BASE[i][1]                        for i in 1:COUNTRY_N]
const COUNTRY_INIT  = [COUNTRY_BASE[i][2]                        for i in 1:COUNTRY_N]
const COUNTRY_RANGE = Tuple{Float64, Float64}[COUNTRY_BASE[i][3] for i in 1:COUNTRY_N]
const COUNTRY_OPTIM = [COUNTRY_BASE[i][5] for i in 1:COUNTRY_N]

# This is used all the time
# Get the modelStart parameter
const COUNTRY_PARAM_START = convert(Int64, findfirst("modelStart" .== COUNTRY_NAMES))
function getCountryStartDate(country::String)::Float64
   return convert(Float64, countryData[country][:params][COUNTRY_PARAM_START])
end

#-------------------------------------------------------------------------------------------------
#--
#-- Define asks indicating variables relating to the disease and to the countries
#-- Note that 2 variables are never optimised to provide an anchor for the others
#-- They are the infectiousness of a symptomatic individual γᵢ and the fatality rate from ICU δᵤ
#--
#PARAMS_diseaseMask = [PARAMS_BASE[i][4]                  for i in 1:PARAMS_N]
#PARAMS_countryMask = ones(length(PARAMS_diseaseMask)) .- PARAMS_diseaseMask

#-------------------------------------------------------------------------------------------------
#--
#-- MODELLING
#--
const COMPARTMENTS = [["S",                 "Susceptible"],
                      ["E",                 "Exposed"],
                      ["I",                 "Infectious"],
                      ["J",                 "Severe"],
                      ["H",                 "Hospitalised"],
                      ["C",                 "ICU"],
                      ["R",                 "Recovered"],
                      ["F",                 "Fully_Recovered"],
                      ["D",                 "Dead"],
                      ["K",                 "Critical"],
                      ["L",                 "Critical_Hospitalised"]]

const COMPARTMENTS_LIST = [v[1] for v in COMPARTMENTS]
const COMPARTMENTS_N    = length(COMPARTMENTS_LIST)

const VARIABLES =  [["BED",                 "Hospital_Beds"],
              ["ICU",                 "ICU_Beds"]]

#--
#-- INDIVIDUAL VARIABLES
#--
const VARIABLES_LIST = [v[1] for v in VARIABLES]
const VARIABLES_N    = length(VARIABLES_LIST)

# Where to find the index of the begiining of the compartment X
function compIndex(C:: String)
   n = findfirst(COMPARTMENTS_LIST .== C) - 1
   return n * nAgeGroup + 1
end

# Where to find the index of the variable X
function varIndex(V:: String)
   n = findfirst(VARIABLES_LIST .== V)
   return COMPARTMENTS_N * nAgeGroup + n
end


text:
__precompile__(true)

module DiffLinearAlgebra

    # Some aliases used repeatedly throughout the package.
    const AV, AM, AVM, AA = AbstractVector, AbstractMatrix, AbstractVecOrMat, AbstractArray
    const SV, SM, SVM, SA = StridedVector, StridedMatrix, StridedVecOrMat, StridedArray
    const AS, ASVM = Union{Real, AA}, Union{Real, AVM}
    const Arg1, Arg2, Arg3 = Type{Val{1}}, Type{Val{2}}, Type{Val{3}}
    const Arg4, Arg5, Arg6 = Type{Val{4}}, Type{Val{5}}, Type{Val{6}}
    const BF = Union{Float32, Float64}
    const DLA = DiffLinearAlgebra

    export ∇, DLA, import_expr

    # Define container and meta-data for defined operations.
    include("util.jl")

    # Define operations, and log them for external use via `ops`.
    include("generic.jl")
    include("blas.jl")
    include("diagonal.jl")
    include("triangular.jl")
    include("uniformscaling.jl")
    include("factorization/cholesky.jl")

end # module


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    FractionPartitioner(fraction, shuffle=true)

A method for partitioning spatial objects according to a given `fraction`.
Optionally `shuffle` elements before partitioning.
"""
struct FractionPartitioner <: AbstractPartitioner
  fraction::Float64
  shuffle::Bool

  function FractionPartitioner(fraction, shuffle)
    @assert 0 < fraction < 1 "fraction must be in interval (0,1)"
    new(fraction, shuffle)
  end
end

FractionPartitioner(fraction) = FractionPartitioner(fraction, true)

function partition(object::AbstractSpatialObject{T,N},
                   p::FractionPartitioner) where {N,T}
  npts = npoints(object)
  frac = round(Int, p.fraction * npts)

  locs = p.shuffle ? randperm(npts) : 1:npts
  subsets = [locs[1:frac], locs[frac+1:npts]]

  SpatialPartition(object, subsets)
end


text:
"""
    inhomogeneous_poisson_process(lam, dt; kwargs...)

Returns a sequence of spike times produced by the inhomogeneous poisson process
with rate `lambda`. `dt` is 1/sampling rate.

To produce a more reliable spike train, increase `sampling_rate_factor`. If
`sampling_rate_factor = k` then the function first computes `kron(lam, ones(k))`
and uses sampling rate `dt/k` (For instance, with the CRCNS data, `k` between 10
and 20 seems to work).

Will not produce spikes with repeated spike times (this could come up if, e.g.
`lam[t]` is very large)

Options:
 * `sampling_rate_factor::Int = 1` Increases sampling rate

Stopping criteria:
 * `Tmax = length(lam) * dt` Will produces spikes until first spike later than this time
 * `max_real_time` Will stop this process early if it's taking forever (i.e. if `lam` is very low). Default value `0` ignores this option.
 * `max_loops` Same as above, counting number of loops instead of amount of real time. Default `0` to ignore.
 * `max_spikes` Same as above, counting number of spikes produced (e.g. if you want the first-so-many spieks produced by a process.)

These additional kwargs can be passed arrays. The first entry will be modified
in place. This allows additional information to be collected about the process
without requiring the function to return multiple values.
 * `total_time` The total run time of the function
 * `total_loops` The total number of times the process looped
 * `total_spikes` The same as `length(times)` where `times` is the output of this function.

This implementation is based on the first algorithm outlined at
http://freakonometrics.hypotheses.org/724

"""
function inhomogeneous_poisson_process(lam::Vector{Float64}, dt::Float64;
    Tmax::Float64=length(lam) * dt, sampling_rate_factor::Int=1,
    max_real_time::Float64=0.0, max_loops::Int=0, max_spikes::Int=0,
    total_time=[], total_loops=[], total_spikes=[],
    exit_status=Symbol[])
    # based on the first algorithm outlined at http://freakonometrics.hypotheses.org/724
    dt = dt/sampling_rate_factor
    # Lam = cumsum_kbn(kron(lam, ones(sampling_rate_factor)) * dt)
    Lam = cumsum_kbn(kron(lam, fill(dt, sampling_rate_factor)))
    times = Vector{Float64}()
    t_last = 0.0
    s = 0.0; t = 0.0
    real_time = 0.0
    n_loops = 0
    n_spikes = 0
    while t < Tmax
        tic()
        u = 1 - rand()
        s = s - log(u)
        idx = searchsortedlast(Lam, s)
        t = idx * dt
        if Tmax > t > t_last
            # Avoid repeated spikes by comparing to t_last
            push!(times, t)
            t_last = t
            n_spikes += 1
        end
        real_time += toq()
        n_loops += 1
        if idx == length(Lam)
            # if we've hit the end of Lam, then we should quit while we're ahead
            push!(exit_status,:end_of_lambda_reached)
            break
        end
        if 0.0 < max_real_time < real_time
            push!(exit_status,:max_real_time_reached)
            break
        end
        if 0 < max_loops < n_loops
            push!(exit_status,:max_loops_reached)
            break
        end
        if 0 < max_spikes < n_spikes
            push!(exit_status,:max_number_of_spikes)
            break
        end
    end
    if t >= Tmax
        push!(exit_status,:max_spike_time_reached)
    end
    if length(total_time) > 0
        total_time[1] = real_time
    end
    if length(total_loops) > 0
        total_loops[1] = n_loops
    end
    if length(total_spikes) > 0
        total_spikes[1] = length(times)
    end
    return times
end

# MAYBEDO maybe implement a version of this function that returns a proper
# SpikeTrains object. Probably not necessary, but keeping this flag here just in
# case.


text:
export svd,
       qr,
       polar,
       eigen,
       factorize

function truncate!(P::Vector{Float64};
                   kwargs...)::Tuple{Float64,Float64}
  maxdim::Int = min(get(kwargs,:maxdim,length(P)), length(P))
  mindim::Int = max(get(kwargs,:mindim,1), 1)
  cutoff::Float64 = max(get(kwargs,:cutoff,0.0), 0.0)
  absoluteCutoff::Bool = get(kwargs,:absoluteCutoff,false)
  doRelCutoff::Bool = get(kwargs,:doRelCutoff,true)

  origm = length(P)
  docut = 0.0

  if P[1]<=0.0
    P[1] = 0.0
    resize!(P,1)
    return 0.,0.
  end

  if origm==1
    docut = P[1]/2
    return 0.,docut
  end

  #Zero out any negative weight
  for n=origm:-1:1
    (P[n] >= 0.0) && break
    P[n] = 0.0
  end

  n = origm
  truncerr = 0.0
  while n > maxdim
    truncerr += P[n]
    n -= 1
  end

  if absoluteCutoff
    #Test if individual prob. weights fall below cutoff
    #rather than using *sum* of discarded weights
    while P[n] <= cutoff && n > mindim
      truncerr += P[n]
      n -= 1
    end
  else
    scale = 1.0
    if doRelCutoff
      scale = sum(P)
      (scale==0.0) && (scale = 1.0)
    end

    #Continue truncating until *sum* of discarded probability
    #weight reaches cutoff reached (or m==mindim)
    while (truncerr+P[n] <= cutoff*scale) && (n > mindim)
      truncerr += P[n]
      n -= 1
    end

    truncerr /= scale
  end

  if n < 1
    n = 1
  end

  if n < origm
    docut = (P[n]+P[n+1])/2
    if abs(P[n]-P[n+1]) < 1E-3*P[n]
      docut += 1E-3*P[n]
    end
  end

  resize!(P,n)

  return truncerr,docut
end

function _permute_for_factorize(A::ITensor,
                                Linds...)
  Ais = inds(A)
  Lis_orig = IndexSet(Linds...)
  Lis = commoninds(Ais,Lis_orig)
  Ris = uniqueinds(Ais,Lis)
  Ais_perm = IndexSet(Lis...,Ris...)
  # TODO: check if hassameinds(Lis,Ais[1:length(Lis)])
  # so that a permute can be avoided
  if inds(A) ≠ Ais_perm
    A = permute(A,Ais_perm)
  end
  return A,Lis,Ris
end

import LinearAlgebra.qr
function qr(A::ITensor,
            Linds...)
  A,Lis,Ris = _permute_for_factorize(A,Linds...)
  Qis,Qstore,Pis,Pstore = storage_qr(store(A),Lis,Ris)
  return ITensor(Qis,Qstore),ITensor(Pis,Pstore)
end

function polar(A::ITensor,
               Linds...)
  A,Lis,Ris = _permute_for_factorize(A,Linds...)
  Qis,Qstore,Pis,Pstore = storage_polar(store(A),Lis,Ris)
  return ITensor(Qis,Qstore),ITensor(Pis,Pstore)
end

import LinearAlgebra.svd
"""
    svd(A::ITensor,
        leftind1::Index,
        leftind2::Index,
        ...
        ;kwargs...)

Singular value decomposition (SVD) of an ITensor A, computed
by treating the "left indices" provided collectively
as a row index, and the remaining "right indices" as a
column index (matricization of a tensor).

Whether the SVD performs a trunction depends on the keyword
arguments provided. The following keyword arguments are recognized:
* `maxdim` [Int]
* `mindim` [Int]
* `cutoff` [Float64]
* `absoluteCutoff` [Bool] Default value: false.
* `doRelCutoff` [Bool] Default value: true.
* `utags` [String] Default value: "Link,u".
* `vtags` [String] Default value: "Link,v".
* `fastSVD` [Bool] Defaut value: false.
"""
function svd(A::ITensor,
             Linds...;
             kwargs...)
  A,Lis,Ris = _permute_for_factorize(A,Linds...)
  Uis,Ustore,Sis,Sstore,Vis,Vstore = storage_svd(store(A),Lis,Ris;kwargs...)

  U = ITensor(Uis,Ustore)
  S = ITensor(Sis,Sstore)
  V = ITensor(Vis,Vstore)
  u = commonindex(U,S)
  v = commonindex(S,V)
  return U,S,V,u,v
end

function _factorize_center(A::ITensor,
                           Linds...;
                           kwargs...)
  tags::TagSet = get(kwargs,:tags,"Link,u")
  U,S,V = svd(A,Linds...;kwargs...)
  u = commonindex(U,S)
  v = commonindex(S,V)
  for ss = 1:dim(u)
    S[ss,ss] = sqrt(S[ss,ss])
  end
  FU = settags(U*S,tags,v)
  FV = settags(S*V,tags,u)
  return FU,FV,commonindex(FU,FV)
end

function _factorize_from_left_svd(A::ITensor,
                                  Linds...;
                                  kwargs...)
  tags::TagSet = get(kwargs,:tags,"Link,u")
  U,S,V = svd(A,Linds...;kwargs...)
  u = commonindex(U,S)
  FU = settags(U,tags,u)
  FV = settags(S*V,tags,u)
  return FU,FV,commonindex(FU,FV)
end

function _factorize_from_right_svd(A::ITensor,
                                   Linds...;
                                   kwargs...)
  tags::TagSet = get(kwargs,:tags,"Link,u")
  U,S,V = svd(A,Linds...;kwargs...)
  v = commonindex(S,V)
  FU = settags(U*S,tags,v)
  FV = settags(V,tags,v)
  return FU,FV,commonindex(FU,FV)
end

function _factorize_from_left_eigen(A::ITensor,
                                    Linds...;
                                    kwargs...)
  A,Lis,Ris = _permute_for_factorize(A,Linds...)
  A² = A*prime(dag(A),Lis)
  FU,D = eigen(A²,Lis,prime(Lis);kwargs...)
  FV = dag(FU)*A
  return FU,FV,commonindex(FU,FV)
end

function _factorize_from_right_eigen(A::ITensor,
                                     Linds...;
                                     kwargs...)
  A,Lis,Ris = _permute_for_factorize(A,Linds...)
  A² = A*prime(dag(A),Ris)
  FV,D = eigen(A²,Ris,prime(Ris); kwargs...)
  FU = A*dag(FV)
  return FU,FV,commonindex(FU,FV)
end

import LinearAlgebra.factorize
"""
factorize(A::ITensor, Linds...; dir = "center", which_factorization = "svd", cutoff = 0.0)

Do a low rank factorization of A either using an SVD or an eigendecomposition of A'A or AA'.
"""
function factorize(A::ITensor,
                   Linds...;
                   kwargs...)
  dir::String = get(kwargs,:dir,"center")
  if dir == "center"
    return _factorize_center(A,Linds...;kwargs...)
  end
  which_factorization::String = get(kwargs,:which_factorization,"svd")
  cutoff::Float64 = get(kwargs,:cutoff,0.0)
  use_eigen = false
  if which_factorization == "eigen" || (which_factorization == "automatic" && cutoff > 1e-12)
    use_eigen = true
  end
  if dir == "fromleft"
    if use_eigen
      return _factorize_from_left_eigen(A,Linds...;kwargs...)
    else
      return _factorize_from_left_svd(A,Linds...;kwargs...)
    end
  elseif dir == "fromright"
    if use_eigen
      return _factorize_from_right_eigen(A,Linds...;kwargs...)
    else
      return _factorize_from_right_svd(A,Linds...;kwargs...)
    end
  end
  error("In factorize, no dir = $dir supported. Use center, fromleft or fromright.")
end

# TODO: add a version that automatically detects the IndexSets
# by matching based on tags
import LinearAlgebra.eigen
function eigen(A::ITensor,
               Linds,
               Rinds;
               kwargs...)
  Lis = IndexSet(Linds)
  Ris = IndexSet(Rinds)
  Ais_perm = IndexSet(Lis...,Ris...)
  !hassameinds(Ais_perm,A) && throw(ErrorException("Input indices must be contained in the ITensor"))
  if inds(A) ≠ Ais_perm
    A = permute(A,Ais_perm)
  end
  #TODO: More of the index analysis should be moved out of storage_eigen
  Uis,Ustore,Dis,Dstore = storage_eigen(store(A),Lis,Ris; kwargs...)
  return ITensor(Uis,Ustore),ITensor(Dis,Dstore)
end


text:
module essentials
	using Gaston
	set(showable="svg")
	using Yao
	using StatsBase: Histogram, fit
	using BitBasis: BitStr
	Θ = LinRange(0, 2π, 100) # 50
	Φ = LinRange(0, π, 20)
	r = 1
	x = [r*cos(θ)*sin(ϕ) for θ in Θ, ϕ in Φ]
	y = [r*sin(θ)*sin(ϕ) for θ in Θ, ϕ in Φ]
	z = [r*cos(ϕ)        for θ in Θ, ϕ in Φ]
	function plot_blochsph(a::Matrix{ComplexF64})
		surf(x, y, z, lc = :turquoise,
			Axes(view = "equal xyz",
				pm3d = "depthorder",
				palette = :gray,
				style = "fill transparent solid 0.1", tics=:false, border=:false, colorbox=:false), w=:pm3d)
		surf!([0,0], [-1,1], [0,0], lc=:red)
		surf!([-1,1], [0,0], [0,0], lc=:brown)
		surf!([0,0], [0,0], [-1,1], lc=:orange)
		surf!([0,0,0,1], [0,0,-1.2,0], [1.2,-1.2,0,0.2], supp=["|0〉","|1〉","x","y"], w="labels")
		p1 = 2 * acos(abs(a[1]))
		p2 = angle(a[2]) - angle(a[1])
		x1 = r*sin(p1)*cos(p2)
		y1 = r*sin(p2)*sin(p1)
		z1 = r*cos(p1)
		surf!([0], [0], [0], supp=[[y1] [-x1] [z1]], w=:vector, lw=5, lc="'turquoise'")
	end

	function plot_hist(x::Array{BitStr{n,Int},1}) where n
		set(preamble="set xtics font ',$(n<=3 ? 15 : 15/(2^(n-3)))'")
		hist = fit(Histogram, Int.(x), 0:2^n)
		bar(hist.edges[1][1:end-1], hist.weights, fc="'dark-red'", Axes(title = :Histogram, xtics = (0:(2^n-1), "|" .* string.(0:(2^n-1), base=2, pad=n) .* "〉"), yrange = "[0:]"))
	end
end

max_stars_repo_name: VarLad/qml-mooc

max_issues_repo_name: VarLad/qml-mooc

max_forks_repo_name: VarLad/qml-mooc


text:
function readimages!(buffer::Matrix{UInt8}, io::IO, index::Integer, nrows::Integer, ncols::Integer)
    seek(io, IMAGEOFFSET + nrows * ncols * (index - 1))
    read!(io, buffer)
end

"""
    readimages(io::IO, index::Integer, nrows::Integer, ncols::Integer)

Jumps to the position of `io` where the bytes for the `index`'th
image are located and reads the next `nrows` * `ncols` bytes. The
read bytes are returned as a `Matrix{UInt8}` of size `(nrows, ncols)`.
"""
function readimages(io::IO, index::Integer, nrows::Integer, ncols::Integer)
    buffer = Array{UInt8}(undef, nrows, ncols)
    readimages!(buffer, io, index, nrows, ncols)
end

"""
    readimages(io::IO, indices::AbstractVector, nrows::Integer, ncols::Integer)

Reads the first `nrows` * `ncols` bytes for each image index in
`indices` and stores them in a `Array{UInt8,3}` of size `(nrows,
ncols, length(indices))` in the same order as denoted by
`indices`.
"""
function readimages(io::IO, indices::AbstractVector, nrows::Integer, ncols::Integer)
    images = Array{UInt8}(undef, nrows, ncols, length(indices))
    buffer = Array{UInt8}(undef, nrows, ncols)
    dst_index = 1
    for src_index in indices
        readimages!(buffer, io, src_index, nrows, ncols)
        copyto!(images, 1 + nrows * ncols * (dst_index - 1), buffer, 1, nrows * ncols)
        dst_index += 1
    end
    images
end

"""
    readimages(file, [indices])

Reads the images denoted by `indices` from `file`. The given
`file` can either be specified using an IO-stream or a string
that denotes the fully qualified path. The conent of `file` is
assumed to be in the MNIST image-file format, as it is described
on the official homepage at http://yann.lecun.com/exdb/mnist/

- if `indices` is an `Integer`, the single image is returned as
  `Matrix{UInt8}` in horizontal major layout, which means that
  the first dimension denotes the pixel *rows* (x), and the
  second dimension denotes the pixel *columns* (y) of the image.

- if `indices` is a `AbstractVector`, the images are returned as
  a 3D array (i.e. a `Array{UInt8,3}`), in which the first
  dimension corresponds to the pixel *rows* (x) of the image, the
  second dimension to the pixel *columns* (y) of the image, and
  the third dimension denotes the index of the image.

- if `indices` is ommited all images are returned
  (as 3D array described above)
"""
function readimages(io::IO, indices)
    _, nimages, nrows, ncols = readimageheader(io)
    @assert minimum(indices) >= 1 && maximum(indices) <= nimages
    readimages(io, indices, nrows, ncols)
end

function readimages(file::AbstractString, index::Integer)
    gzopen(file, "r") do io
        readimages(io, index)
    end::Matrix{UInt8}
end

function readimages(file::AbstractString, indices::AbstractVector)
    gzopen(file, "r") do io
        readimages(io, indices)
    end::Array{UInt8,3}
end

function readimages(file::AbstractString)
    gzopen(file, "r") do io
        _, nimages, nrows, ncols = readimageheader(io)
        readimages(io, 1:nimages, nrows, ncols)
    end::Array{UInt8,3}
end