text:
module TestSelect

using FactCheck
using DataCubes

immutable Wrap{T}
  elem::T
end

facts("Select tests") do
  col1 = nalift(reshape(500:-1:1, 10, 50))
  col2 = nalift(reshape(1.0*(1:500), 10, 50))
  col3 = nalift(reshape(map(i->string("sym_",i), 1:500), 10, 50))
  col4 = nalift(hcat(rand(10,30), fill(:testsym, 10, 20)))
  axis1c1 = DictArray(k1=nalift(collect(101:110)), k2=nalift(collect(201:210)))
  axis1c2 = DictArray(r1=nalift(map(i->string("a_",i), 1:50)))
  d = DictArray(c1=col1, c2=col2, c3=col3, c4=col4)
  lar = LabeledArray(d, axis1=axis1c1, axis2=axis1c2)
  context("@select tests") do
    @fact @select(lar) === lar --> true
    context("aggregate tests") do
      @fact @select(lar,:c2,:c3).data --> pick(lar,[:c2,:c3])
      @fact @select(lar,c2=broadcast(*, _c2, 2),:c3).data --> @darr(c2=2*pick(lar,:c2),c3=pick(lar,:c3))
      @fact @select(lar,c2=broadcast(*, _c2, 2),:c3) --> @select(lar,c2=broadcast(*, _[:c2], 2),:c3)
      @fact @select(larr(a=[1,2,3],b=[4,5,6]), :a=>sum(_a), :b=>length(_b)) --> nalift(LDict(:a=>6,:b=>3))
      @fact @select(@darr(a=[1 2 3 NA]), s=sum(_a), where[broadcast(>=, _a,2.0)]) --> nalift(LDict(:s=>5))
      @fact @select(larr(a=[1 2 3;4 5 6]), :a) --> larr(a=[1 2 3;4 5 6])
      @fact_throws @select(larr(rand(3,4)), :a)
      @fact_throws @select(larr(rand(5,3)), a=_a)
    end
    context("condition tests") do
      @fact @select(lar, where[broadcast(<, 10, _c2) & broadcast(<, _c2, 25)]) --> begin
        base = @larr(c1=col1[1:10,2:3],c2=col2[1:10,2:3],c3=col3[1:10,2:3],c4=col4[1:10,2:3],axis1[k1=collect(101:110),k2=collect(201:210)],axis2[r1=["a_2","a_3"]])
        dcube.setna!(base, 5:10, 2)
        base
      end
      @fact @select(lar, c1=broadcast(*, _c1, 2), where[broadcast(<, 10, _c2) & broadcast(<, _c2, 25)]) --> begin
        base = @larr(c1=broadcast(*, 2, col1[1:10,2:3]),axis1[k1=collect(101:110),k2=collect(201:210)],axis2[r1=["a_2","a_3"]])
        dcube.setna!(base, 5:10, 2)
        base
      end
    end
    context("by tests") do
      @fact size(@select(lar, by[k1=_k2 - _c1,:r1], by[m1=broadcast(>, _k2, 205)])) --> (500,2)
      @fact size(@select(lar, by[k1=_k2 - _k1], by[m1=broadcast(>, _k2, 205)])) --> (1,2)
      @fact @select(lar, by[k1=_k2 - _k1], by[m1=broadcast(>, _k2, 205)], c1=sum(igna(_c1))) --> @larr(c1=[63250 62000],axis1[k1=[100]],axis2[m1=[false,true]])
      f1 = fld -> sum(igna(fld))
      @fact @select(lar, by[k1=_k2 - _k1], by[m1=broadcast(>, _k2, 205)], c1=f1(_c1)) --> @larr(c1=[63250 62000],axis1[k1=[100]],axis2[m1=[false,true]])
      f2 = d -> sum(igna(d[:c1]))
      @fact @select(lar, by[k1=_k2 - _k1], by[m1=broadcast(>, _k2, 205)], c1=f2(_)) --> @larr(c1=[63250 62000],axis1[k1=[100]],axis2[m1=[false,true]])
      @fact size(@select(lar, by[k1=_k2 - _c1,:r1], by[m1=broadcast(>, _k2, 205)], c1= broadcast(*, _c1, _c2), :c2)) --> (500,2)
      @fact @select(@larr(a=collect(1:100),b=enumeration(repmat(collect(1:10),10))), by[:b], :a=>sum(_a), where[broadcast(>, _a, 56)]) --> sort(@larr(a=[385,390,395,400,304,308,312,316,320,324],axis1[b=[7,8,9,10,1,2,3,4,5,6]]), 1, :b)
      @fact @select(@larr(a=collect(1:100),b=repmat(collect(1:10),10), c=repmat(collect(1:10),10)), by[:b,:c], :a=>sum(_a), where[broadcast(>, _a, 56)]) --> sort(@larr(a=[385,390,395,400,304,308,312,316,320,324],axis1[b=[7,8,9,10,1,2,3,4,5,6],c=[7,8,9,10,1,2,3,4,5,6]]), 1, :b)
      @fact @select(@larr(a=collect(1:100),b=enumeration(repmat(collect(1:10),10)), c=repmat(collect(1:10),10)), by[:b,:c], :a=>sum(_a), where[broadcast(>, _a, 56)]) --> sort(@larr(a=[385,390,395,400,304,308,312,316,320,324],axis1[b=[7,8,9,10,1,2,3,4,5,6],c=[7,8,9,10,1,2,3,4,5,6]]), 1, :b)
      @fact @select(@larr(a=collect(1:100),b=repmat(collect(1:10),10), c=enumeration(repmat(collect(1:10),10))), by[:b,:c], :a=>sum(_a), where[broadcast(>, _a, 56)]) --> sort(@larr(a=[385,390,395,400,304,308,312,316,320,324],axis1[b=[7,8,9,10,1,2,3,4,5,6],c=[7,8,9,10,1,2,3,4,5,6]]), 1, :b)
      @fact @select(@larr(a=collect(1:100),b=enumeration(repmat(collect(1:10),10)), c=enumeration(repmat(collect(1:10),10))), by[:b,:c], :a=>sum(_a), where[broadcast(>, _a, 56)]) --> sort(@larr(a=[385,390,395,400,304,308,312,316,320,324],axis1[b=[7,8,9,10,1,2,3,4,5,6],c=[7,8,9,10,1,2,3,4,5,6]]), 1, :b)
    end
  end
  context("selct tests") do
    context("aggregate tests") do
      @fact selct(lar,:c2,:c3).data --> pick(lar,[:c2,:c3])
      @fact selct(lar,c2=d->broadcast(*, d[:c2], 2),:c3).data --> reorder(@darr(c2=2*pick(lar,:c2),c3=pick(lar,(:c3,))[1]), :c3)
      @fact selct(lar,c2=d->broadcast(*, d[:c2], 2),:c3) --> reorder(@select(lar,c2=broadcast(*, _[:c2], 2),:c3), :c3)
    end
    context("condition tests") do
      @fact selct(lar, where=d->broadcast(<, 10.0, d[:c2]) & broadcast(<, d[:c2], 25)) --> begin
        base = @larr(c1=col1[1:10,2:3],c2=col2[1:10,2:3],c3=col3[1:10,2:3],c4=col4[1:10,2:3],axis1[k1=collect(101:110),k2=collect(201:210)],axis2[r1=["a_2","a_3"]])
        dcube.setna!(base, 5:10, 2)
        base
      end
      @fact selct(lar, c1=d->broadcast(*, d[:c1], 2), where=Any[d->broadcast(<, 10.0, d[:c2]) & broadcast(<, d[:c2], 25)]) --> begin
        base = @larr(c1=broadcast(*, 2, col1[1:10,2:3]),axis1[k1=collect(101:110),k2=collect(201:210)],axis2[r1=["a_2","a_3"]])
        dcube.setna!(base, 5:10, 2)
        base
      end
    end
    context("by tests") do
      @fact size(selct(lar, by=Any[:k1=>d->d[:k2]-d[:c1],:r1], by=:m1=>d->broadcast(>, d[:k2], 205))) --> (500,2)
      @fact size(selct(lar, by=:k1=>d->d[:k2]-d[:k1], by=Any[:m1=>d->broadcast(>, d[:k2], 205)])) --> (1,2)
      @fact selct(lar, by=:k1=>d->d[:k2] - d[:k1], by=Any[:m1=>d->broadcast(>, d[:k2], 205)], c1=d->sum(igna(d[:c1]))) --> @larr(c1=[63250 62000],axis1[k1=[100]],axis2[m1=[false,true]])
      f1 = fld -> sum(igna(fld))
      @fact selct(lar, by=:k1=>d->d[:k2]-d[:k1], by=:m1=>d->broadcast(>, d[:k2], 205), c1=d->f1(d[:c1])) --> @larr(c1=[63250 62000],axis1[k1=[100]],axis2[m1=[false,true]])
      f2 = d -> sum(igna(d[:c1]))
      @fact selct(lar, by=Any[:k1=>d->d[:k2]-d[:k1]], by=:m1=>d->broadcast(>, d[:k2], 205), c1=d->f2(d)) --> @larr(c1=[63250 62000],axis1[k1=[100]],axis2[m1=[false,true]])
      @fact size(selct(lar, by=Any[:k1=>d->d[:k2]-d[:c1],:r1], by=:m1=>d->broadcast(>, d[:k2], 205), :c1=>d->broadcast(*, d[:c1], d[:c2]), :c2)) --> (500,2)
      @fact selct(@larr(a=collect(1:100),b=repmat(collect(1:10),10)), by=[:b], :a=>d->sum(d[:a]), where=d->broadcast(>, d[:a], 56)) --> sort(@larr(a=[385,390,395,400,304,308,312,316,320,324],axis1[b=[7,8,9,10,1,2,3,4,5,6]]), 1, :b)
      @fact selct(@darr(a=collect(1:100),b=repmat(collect(1:10),10)), by=[:b], :a=>d->sum(d[:a]), where=d->broadcast(>, d[:a], 56)) --> sort(@larr(a=[385,390,395,400,304,308,312,316,320,324],axis1[b=[7,8,9,10,1,2,3,4,5,6]]), 1, :b)
    end
  end

  context("@update tests") do
    context("aggregate tests") do
      @fact @update(lar,c2=broadcast(*, _c2, 2)) --> (temp=copy(lar);temp.data.data[:c2]=broadcast(*, temp.data.data[:c2], 2);temp)
      @fact @update(lar,c2=broadcast(*, _[:c2], 2)) --> (temp=copy(lar);temp.data.data[:c2]=broadcast(*, temp.data.data[:c2], 2);temp)
    end
    context("condition tests") do
      @fact @update(lar, c1=broadcast(*, _c1, 2), where[broadcast(<, 10, _c2) & broadcast(<, _c2, 25)]) --> begin
        temp = copy(lar)
        inds = ignabool(broadcast(<, 10, temp.data.data[:c2]) & broadcast(<, temp.data.data[:c2], 25))
        temp.data.data[:c1][inds] = broadcast(*, 2, temp.data.data[:c1][inds])
        temp
      end
    end
  end
  context("update tests") do
    context("aggregate tests") do
      @fact update(lar,c2=d->broadcast(*, d[:c2], 2)) --> (temp=copy(lar);temp.data.data[:c2]=broadcast(*, temp.data.data[:c2], 2);temp)
    end
    context("condition tests") do
      @fact update(lar, c1=d->broadcast(*, d[:c1], 2), where=d->broadcast(<,10,d[:c2]) & broadcast(<,d[:c2], 25)) --> begin
        temp = copy(lar)
        inds = ignabool(broadcast(<, 10, temp.data.data[:c2]) & broadcast(<, temp.data.data[:c2], 25))
        temp.data.data[:c1][inds] = broadcast(*, 2, temp.data.data[:c1][inds])
        temp
      end
    end
  end

  context("conversion tests") do
    @fact @select(@darr(a=[1,2,3]),b=[4,5,6]) --> @darr(b=[4,5,6])
    @fact @update(@darr(a=[1,2,3]),b=[4,5,6]) --> @darr(a=[1,2,3],b=[4,5,6])
    @fact @select(@larr(a=[1,2,3]),b=[4,5,6]) --> @larr(b=[4,5,6])
    @fact @update(@larr(a=[1,2,3]),b=[4,5,6]) --> @larr(a=[1,2,3],b=[4,5,6])
  end
  context("@select misc tests") do
    @fact @select(@larr(a=[1 2 3;4 5 6],b=[10 11 12;13 14 15],c=[:a :b :b;:b :a :a]), k=_[:a], where[broadcast(>, _a, 2)]) --> @larr(k=[NA NA 3;4 5 6])
    @fact @select(@larr(a=[1 2 3;4 5 6],b=[10 11 12;13 14 15],c=[:a :b :b;:b :a :a]), k=_[_c]) -->@larr(k=[1 11 12;13 5 6])
    @fact @select(@larr(a=[1 2 3;4 5 6],b=[10 11 12;13 14 15],c=[:a :b :b;:b :a :a]), k=_[_c], where[broadcast(>, _a, 2)]) --> @larr(k=[NA NA 12;13 5 6])
    # test whether sortable byvariables are sorted.
    @fact @select(@larr(a=[1 2 3 4 5],b=[:a :b :a :b :b]), by[:b], r=sum(_a)) --> @larr(axis1[b=[:a,:b]], r=[4,11])
    @fact @select(@larr(a=[1 2 3 4 5],b=[:b :a :a :b :b]), by[:b], r=sum(_a)) --> @larr(axis1[b=[:a,:b]], r=[5,10])
    # test whether unsortable byvariables are unsorted.
    @fact @select(@larr(a=[1 2 3 4 5],b=[Wrap(:a) Wrap(:b) Wrap(:a) Wrap(:b) Wrap(:b)]), by[:b], r=sum(_a)) --> @larr(axis1[b=[Wrap(:a),Wrap(:b)]], r=[4,11])
    @fact @select(@larr(a=[1 2 3 4 5],b=[Wrap(:b) Wrap(:a) Wrap(:a) Wrap(:b) Wrap(:b)]), by[:b], r=sum(_a)) --> @larr(axis1[b=[Wrap(:b),Wrap(:a)]], r=[10,5])
    @fact size(@select(@larr(a=fill(:sym,10),b=[1,1,1,2,2,3,3,3,3,3],c=[1,1,2,2,2,2,2,3,3,3]), by[:b],by[:c])) --> (3,3)
    @fact isna(@select(@larr(a=fill(:sym,10),b=[1,1,1,2,2,3,3,3,3,3],c=[1,1,2,2,2,2,2,3,3,3]), by[:b],by[:c])) --> [false false true;true false true;true false false]
    @fact @select(@larr(a=fill(:sym,10),b=[1,1,1,2,2,3,3,3,3,3],c=[1,1,2,2,2,2,2,3,3,3]), by[:b],by[:c], a=length(_a)) --> @larr(a=[2 1 NA;NA 2 NA;NA 2 3], axis1[b=[1,2,3]], axis2[c=[1,2,3]])
    @fact @select(larr(a=[1 2 3;4 5 6], b=['a' 'b' 'c';'d' 'e' 'f']), c=length(_), d=_a) --> reshape(larr(c=fill(6,6), d=[1,4,2,5,3,6]), 2, 3)
    @fact @select(larr(a=[1 2 3;4 5 6], b=['a' 'b' 'c';'d' 'e' 'f']), c=length(_), d=_a, where[broadcast(>, _a, 2)]) --> reshape(@larr(c=[NA,4,NA,4,4,4], d=[NA,4,NA,5,3,6]), 2, 3)
    @fact @select(larr([1 2;3 4;5 6],axis1=darr(k=[10,11,12])), where[broadcast(==, _k, 11)]) --> larr([3 4], axis1=darr(k=[11]))
  end
  context("selct misc tests") do
    @fact selct(@larr(a=[1 2 3;4 5 6],b=[10 11 12;13 14 15],c=[:a :b :b;:b :a :a]), k=d->d[:a], where=d->broadcast(>, d[:a], 2)) --> @larr(k=[NA NA 3;4 5 6])
    @fact selct(@larr(a=[1 2 3;4 5 6],b=[10 11 12;13 14 15],c=[:a :b :b;:b :a :a]), k=d->d[:a], where=Any[d->broadcast(>, d[:a], 2)]) --> @larr(k=[NA NA 3;4 5 6])
    @fact selct(@larr(a=[1 2 3;4 5 6],b=[10 11 12;13 14 15],c=[:a :b :b;:b :a :a]), k=d->d[:a], where=Any[d->broadcast(>, d[:a], 2)]) --> @larr(k=[NA NA 3;4 5 6])
    @fact selct(@larr(a=[1 2 3;4 5 6],b=[10 11 12;13 14 15],c=[:a :b :b;:b :a :a]), k=d->d[d[:c]]) -->@larr(k=[1 11 12;13 5 6])
    @fact selct(@larr(a=[1 2 3;4 5 6],b=[10 11 12;13 14 15],c=[:a :b :b;:b :a :a]), k=d->d[d[:c]], where=d->broadcast(>, d[:a], 2)) --> @larr(k=[NA NA 12;13 5 6])
    # test whether sortable byvariables are sorted.
    @fact selct(@larr(a=[1 2 3 4 5],b=[:a :b :a :b :b]), by=[:b], r=d->sum(d[:a])) --> @larr(axis1[b=[:a,:b]], r=[4,11])
    @fact selct(@larr(a=[1 2 3 4 5],b=[:b :a :a :b :b]), by=:b, r=d->sum(d[:a])) --> @larr(axis1[b=[:a,:b]], r=[5,10])
    # test whether unsortable byvariables are unsorted.
    @fact selct(@larr(a=[1 2 3 4 5],b=[Wrap(:a) Wrap(:b) Wrap(:a) Wrap(:b) Wrap(:b)]), by=[:b], r=d->sum(d[:a])) --> @larr(axis1[b=[Wrap(:a),Wrap(:b)]], r=[4,11])
    @fact selct(@larr(a=[1 2 3 4 5],b=[Wrap(:b) Wrap(:a) Wrap(:a) Wrap(:b) Wrap(:b)]), by=:b, r=d->sum(d[:a])) --> @larr(axis1[b=[Wrap(:b),Wrap(:a)]], r=[10,5])
    @fact size(selct(@larr(a=fill(:sym,10),b=[1,1,1,2,2,3,3,3,3,3],c=[1,1,2,2,2,2,2,3,3,3]), by=:b,by=:c)) --> (3,3)
    @fact isna(selct(@larr(a=fill(:sym,10),b=[1,1,1,2,2,3,3,3,3,3],c=[1,1,2,2,2,2,2,3,3,3]), by=:b,by=[:c])) --> [false false true;true false true;true false false]
    @fact selct(@larr(a=fill(:sym,10),b=[1,1,1,2,2,3,3,3,3,3],c=[1,1,2,2,2,2,2,3,3,3]), by=[:b],by=:c, :a=>d->length(d[:a])) --> @larr(a=[2 1 NA;NA 2 NA;NA 2 3], axis1[b=[1,2,3]], axis2[c=[1,2,3]])
    @fact selct(larr(a=[1 2 3;4 5 6], b=['a' 'b' 'c';'d' 'e' 'f']), :c=>d->length(d), d=d->d[:a]) --> reshape(larr(c=fill(6,6), d=[1,4,2,5,3,6]), 2, 3)
    # note the reversal of order because the keyword arguments come after pair arguments.
    @fact selct(larr(a=[1 2 3;4 5 6], b=['a' 'b' 'c';'d' 'e' 'f']), c=d->length(d), :d=>d->d[:a], where=d->broadcast(>, d[:a], 2)) --> reshape(@larr(d=[NA,4,NA,5,3,6],c=[NA,4,NA,4,4,4]), 2, 3)
    @fact selct(larr([1 2;3 4;5 6],axis1=darr(k=[10,11,12])), where=Any[d->broadcast(==, d[:k], 11)]) --> larr([3 4], axis1=darr(k=[11]))
  end
  context("misc @update misc tests") do
    @fact @update(larr(a=collect(1.0*(1:10)), b=[1,1,1,2,2,3,3,3,3,3]),a=mean(_a), by[:b]) --> larr(a=[2.0,2.0,2.0,4.5,4.5,8.0,8.0,8.0,8.0,8.0], b=[1,1,1,2,2,3,3,3,3,3])
    # TODO a=1:10 was fine in julia v0.5. In v0.6, the int type complains if we update it using a = mean(_a).
    @fact @update(larr(a=1.0*1:10, b=[1,1,1,2,2,3,3,3,3,3]),a=mean(_a), by[:b]) --> larr(a=[2.0,2.0,2.0,4.5,4.5,8.0,8.0,8.0,8.0,8.0], b=[1,1,1,2,2,3,3,3,3,3])
    @fact @update(larr(a=collect(1.0*(1:10)), b=[1,1,1,2,2,3,3,3,3,3], c=[11,12,13,14,11,12,13,14,11,12]),a=mean(_a), by[:b,:c]) --> larr(a=[1.0,2.0,3.0,4.0,5.0,8.0,7.0,8.0,9.0,8.0], b=[1,1,1,2,2,3,3,3,3,3], c=[11,12,13,14,11,12,13,14,11,12])
    @fact @update(larr(a=1.0*1:10, b=[1,1,1,2,2,3,3,3,3,3], c=[11,12,13,14,11,12,13,14,11,12]),a=mean(_a), by[:b,:c]) --> larr(a=[1.0,2.0,3.0,4.0,5.0,8.0,7.0,8.0,9.0,8.0], b=[1,1,1,2,2,3,3,3,3,3], c=[11,12,13,14,11,12,13,14,11,12])
    @fact @update(larr(a=collect(1.0*(1:10)), b=[1,1,1,2,2,3,3,3,3,3], c=[11,12,13,14,11,12,13,14,11,12]),a=mean(_a), by[:b,:c]) --> @update(larr(a=collect(1.0*(1:10)), b=[1,1,1,2,2,3,3,3,3,3], c=[11,12,13,14,11,12,13,14,11,12]),a=mean(_a), by[:b], by[:c])
    @fact @update(larr(a=collect(1.0*(1:10)), b=[1,1,1,2,2,3,3,3,3,3], c=[11,12,13,14,11,12,13,14,11,12]),ma=mean(_a), by[:b,:c]) --> larr(a=collect(1.0*(1:10)), b=[1,1,1,2,2,3,3,3,3,3], c=[11,12,13,14,11,12,13,14,11,12], ma=1.0*[1,2,3,4,5,8,7,8,9,8])
    @fact @update(larr(a=collect(1.0*(1:10)), b=[1,1,1,2,2,3,3,3,3,3], c=[11,12,13,14,11,12,13,14,11,12]),ma=reverse(_a), by[:b,:c]) --> larr(a=collect(1.0*(1:10)), b=[1,1,1,2,2,3,3,3,3,3], c=[11,12,13,14,11,12,13,14,11,12], ma=1.0*[1,2,3,4,5,10,7,8,9,6])
    @fact @update(larr(a=collect(1.0*(1:10)), b=[1,1,1,2,2,3,3,3,3,3],c=[11,12,13,14,11,12,13,14,11,12]),a=mean(_a), by[:b,:c], where[broadcast(<=, _c, 12)]) --> larr(a=1.0*[1,2,3,4,5,8,7,8,9,8],b=[1,1,1,2,2,3,3,3,3,3],c=[11,12,13,14,11,12,13,14,11,12])
    @fact @update(larr(a=collect(1.0*(1:10)),c=[11,12,13,14,11,12,13,14,11,12]),a=reverse(_a), by[:c], where[broadcast(<=, _c, 12)]) --> larr(a=1.0*[9,10,3,4,5,6,7,8,1,2], c=[11,12,13,14,11,12,13,14,11,12])
    @fact @update(larr(a=1.0*1:10,c=enumeration([11,12,13,14,11,12,13,14,11,12])),a=1.0*reverse(_a), by[:c], where[broadcast(<=, _c, 12)]) --> larr(a=1.0*[9,10,3,4,5,6,7,8,1,2], c=[11,12,13,14,11,12,13,14,11,12])
    @fact @update(larr(a=collect(1.0*(1:10)),c=enumeration([11,12,13,14,11,12,13,14,11,12])),ma=mean(_a), a=reverse(_a), by[:c], where[broadcast(<=, _c, 12)]) --> @larr(a=1.0*[9,10,3,4,5,6,7,8,1,2], c=[11,12,13,14,11,12,13,14,11,12], ma=[5.0,6.0,NA,NA,5.0,6.0,NA,NA,5.0,6.0])
    @fact @update(darr(a=1.0*[1,2,3,4,5],b=[1,1,2,2,2]),by[b=broadcast(*, _b, 2)],a=reverse(_a),where[broadcast(<, _a, 5)]) --> darr(a=[2,1,4,3,5],b=[1,1,2,2,2])
    @fact @update(darr(a=1.0*[1,2,3,4,5],b=enumeration([1,1,2,2,2])),by[b=broadcast(*, _b, 2)],a=mean(_a),where[broadcast(<, _a, 5)]) --> darr(a=[1.5,1.5,3.5,3.5,5.0],b=[1,1,2,2,2])
    @fact @update(darr(a=[1,2,3,4,5],b=[1,1,2,2,2]),by[b=broadcast(*, _b, 2)],c=mean(_a),where[broadcast(<, _a, 5)]) --> reorder(@darr(c=[1.5,1.5,3.5,3.5,NA],b=[1,1,2,2,2],a=[1,2,3,4,5]), :a,:b,:c)
    @fact @update(darr(a=[1,2,3,4,5],b=[1,1,2,2,2]),c=_a*2.0,where[broadcast(<, _a, 5)]) --> darr(a=[1,2,3,4,5],b=[1,1,2,2,2],c=@nalift([2.0,4.0,6.0,8.0,NA]))
    @fact @update(darr(a=[1,2,3,4,5],b=[1,1,2,2,2]),c=mean(_a),where[broadcast(<, _a, 5)]) --> darr(a=[1,2,3,4,5],b=[1,1,2,2,2],c=@nalift([2.5,2.5,2.5,2.5,NA]))
    # Let's think about this is okay, where we just return the original labeled array when there is nothing to update.
    @fact @update(darr(a=[1,2,3,4,5],b=[1,1,2,2,2]),c=mean(_a),where[broadcast(<, _a, 5)],where[broadcast(>, _a, 5)]) --> darr(a=[1,2,3,4,5],b=[1,1,2,2,2])
    @fact @update(darr(a=[1,2,3,4,5],b=[1,1,2,2,2])) --> darr(a=[1,2,3,4,5],b=[1,1,2,2,2])
    @fact update(@larr(a=1:10, b=[1,2,3,NA,NA,NA,1,1,2,3], c=[:x,:x,:x,:x,:y,:y,:y,:z,:z,:z]), a=d->sum(d[:a]), d=d->reverse(broadcast(*, d[:a], d[:b])), where=[d-> ~isna(d[:b])], by=[:b]) --> @larr(a=[16,11,13,4,5,6,16,16,11,13],b=[1,2,3,NA,NA,NA,1,1,2,3],c=[:x,:x,:x,:x,:y,:y,:y,:z,:z,:z],d=[8,18,30,NA,NA,NA,7,1,4,9])
    @fact @select(larr(a=[1,1,2,3,4],b=enumeration([:a,:a,:b,:b,:b])), ct=length(_), by[:b]) --> larr(axis1=darr(b=[:a,:b]), ct=[2,3])
    @fact typeof(pickaxis(@select(larr(a=[1,1,2,3,4],b=enumeration([:a,:a,:b,:b,:b])), ct=length(_), by[:b]),1)) --> DataCubes.DictArray{Symbol,1,DataCubes.AbstractArrayWrapper{Nullable{Symbol},1,DataCubes.EnumerationArray{Symbol,1,Array{Int64,1},Int64}},Nullable{Symbol}}
    @fact @select(larr(a=[1,1,2,3,4],b=enumeration([1.0,1.0,2.0,1.0,2.0])), ct=length(_), by[:b]) --> @larr(axis1[b=[1.0,2.0]], ct=[3,2])
    @fact (@rap typeof pickaxis(_,1) @select(larr(a=[1,1,2,3,4],b=enumeration([1.0,1.0,2.0,1.0,2.0])), ct=length(_), by[:b])) --> DataCubes.DictArray{Symbol,1,DataCubes.AbstractArrayWrapper{Nullable{Float64},1,DataCubes.FloatNAArray{Float64,1,Array{Float64,1}}},Nullable{Float64}}
    @fact @select(larr(a=[1 2 3;4 5 6], axis2=darr(r=enumeration([:a,:a,:b]))), by[:r], s=sum(_a)) --> @larr(axis1[r=[:a,:b]], s=[12,9])
    @fact typeof(pickaxis(@select(larr(a=[1 2 3;4 5 6], axis2=darr(r=enumeration([:a,:a,:b]))), by[:r], s=sum(_a)),1)) --> DataCubes.DictArray{Symbol,1,DataCubes.AbstractArrayWrapper{Nullable{Symbol},1,DataCubes.EnumerationArray{Symbol,1,Array{Int64,1},Int64}},Nullable{Symbol}}
    @fact @select(larr(a=[1 2 3;4 5 6], axis2=darr(r=enumeration([:a,:a,:b]))), by[r=map(x->string(x.value),_r)], s=sum(_a)) --> @larr(axis1[r=["a","b"]], s=[12,9])
    @fact typeof(pickaxis(@select(larr(a=[1 2 3;4 5 6], axis2=darr(r=enumeration([:a,:a,:b]))), by[r=map(x->string(x.value),_r)], s=sum(_a)),1)) --> DataCubes.DictArray{Symbol,1,DataCubes.AbstractArrayWrapper{Nullable{String},1,Array{Nullable{String},1}},Nullable{String}}
    @fact @select(larr(a=[1 2 3;4 5 6], axis2=darr(r=[1.0,1.0,2.0])), by[:r], s=sum(_a)) --> @larr(axis1[r=[1.0,2.0]], s=[12,9])
    @fact typeof(pickaxis(@select(larr(a=[1 2 3;4 5 6], axis2=darr(r=[1.0,1.0,2.0])), by[:r], s=sum(_a)),1)) --> DataCubes.DictArray{Symbol,1,DataCubes.AbstractArrayWrapper{Nullable{Float64},1,DataCubes.FloatNAArray{Float64,1,Array{Float64,1}}},Nullable{Float64}}
    @fact @select(larr(a=[1 2 3;4 5 6], axis2=darr(r=[1.0,1.0,2.0])), by[r=broadcast(*, _r, 2)], s=sum(_a)) --> @larr(axis1[r=[2.0,4.0]], s=[12,9])
    @fact typeof(pickaxis(@select(larr(a=[1 2 3;4 5 6], axis2=darr(r=[1.0,1.0,2.0])), by[r=broadcast(*,_r, 2)], s=sum(_a)),1)) --> DataCubes.DictArray{Symbol,1,DataCubes.AbstractArrayWrapper{Nullable{Float64},1,DataCubes.FloatNAArray{Float64,1,Array{Float64,1}}},Nullable{Float64}}
    @fact typeof(view(@select(larr(a=rand(5,3),axis2=darr(r=enumeration([:X,:Y,:Z]))),:r),2,1:2).data.data[:r]) -->  typeof(view(@select(larr(a=rand(5,3),axis2=darr(r=enumeration([:X,:Y,:Z]))),:r),2,1:2).data.data[:r])
    @fact @select(@larr(a=[1 2 3;4 5 6], axis1[k=[:a,:b]], axis2[r=['X','Y','Z']]), by[:r], s=sum(_a)) --> larr(axis1=darr(r=['X','Y','Z']), s=[5,7,9])
    @fact @select(@larr(a=[1 2 3;4 5 6], axis1[k=[:a,:b]], axis2[r=['X','Y','Z']]), by[:r,:r1=>_r], s=sum(_a)) --> larr(axis1=darr(r=['X','Y','Z'],r1=['X','Y','Z']), s=[5,7,9])
    @fact @select(@larr(a=[1 2 3;4 5 6], axis1[k=[:a,:b]], axis2[r=['X','Y','Z']]), by[:r,:r1=>_r,r2=_r], s=sum(_a)) --> larr(axis1=darr(r=['X','Y','Z'],r1=['X','Y','Z'],r2=['X','Y','Z']), s=[5,7,9])
    @fact @select(@larr(a=[1 2 3;4 5 6], axis1[k=[:a,:b]], axis2[r=['X','Y','Z']]), by[:r,:r1=>_r,r2=_r,:r3=>_r], s=sum(_a)) --> larr(axis1=darr(r=['X','Y','Z'],r1=['X','Y','Z'],r2=['X','Y','Z'],r3=['X','Y','Z']), s=[5,7,9])
    @fact @select(@larr(a=[1 2 3;4 5 6], axis1[k=[:a,:b]], axis2[r=['X','Y','Z']]), by[:r,:r1=>_r,r2=_r,:r3=>_r,:r4=>_r], s=sum(_a)) --> larr(axis1=darr(r=['X','Y','Z'],r1=['X','Y','Z'],r2=['X','Y','Z'],r3=['X','Y','Z'],r4=['X','Y','Z']), s=[5,7,9])
    @fact @select(@larr(a=[1 2 3;4 5 6], axis1[k=[:a,:b]], axis2[r=['X','Y','Z']]), by[:r,:r1=>_r,r2=_r,:r3=>_r,:r4=>_r,r5=_r], s=sum(_a)) --> larr(axis1=darr(r=['X','Y','Z'],r1=['X','Y','Z'],r2=['X','Y','Z'],r3=['X','Y','Z'],r4=['X','Y','Z'],r5=['X','Y','Z']), s=[5,7,9])
    @fact @select(@larr(a=[1 2 3;4 5 6], axis1[k=[:a,:b]], axis2[r=['X','Y','Z']]), by[:r,:r1=>_r,r2=_r,:r3=>_r,:r4=>_r,r5=_r,r6=_r], s=sum(_a)) --> larr(axis1=darr(r=['X','Y','Z'],r1=['X','Y','Z'],r2=['X','Y','Z'],r3=['X','Y','Z'],r4=['X','Y','Z'],r5=['X','Y','Z'],r6=['X','Y','Z']), s=[5,7,9])
    @fact @select(@larr(a=[1 2 3;4 5 6], axis1[k=[:a,:b]], axis2[r=['X','Y','Z']]), by[:r,:r1=>_r,r2=_r,:r3=>_r,:r4=>_r,r5=_r,r6=_r,r7=_r], s=sum(_a)) --> larr(axis1=darr(r=['X','Y','Z'],r1=['X','Y','Z'],r2=['X','Y','Z'],r3=['X','Y','Z'],r4=['X','Y','Z'],r5=['X','Y','Z'],r6=['X','Y','Z'],r7=['X','Y','Z']), s=[5,7,9])
    @fact size(@select(darr(a=repmat(1:100,100),b=1.0*1:10_000), by[:a], s=msum(_b))) --> (100,)
    @fact size(@update(darr(a=repmat(1:100,100),b=1.0*1:10_000), by[:a], s=msum(_b))) --> (10000,)
    @fact size(@select(darr(a=repmat(1:100,100),b=1.0*1:10000), by[:a], s=mprod(_b))) --> (100,)
    @fact size(@update(darr(a=repmat(1:100,100),b=1.0*1:10000), by[:a], s=mprod(_b, rev=true, window=10))) --> (10000,)
    @fact size(@select(darr(a=repmat(1:100,100),b=1.0*1:10000), by[:a], s=mmean(_b, window=10))) --> (100,)
    @fact size(@update(darr(a=repmat(1:100,100),b=1.0*1:10000), by[:a], s=mmean(_b, window=1))) --> (10000,)
    @fact size(@select(darr(a=repmat(1:100,100),b=1.0*1:10000), by[:a], s=mminimum(_b))) --> (100,)
    @fact size(@update(darr(a=repmat(1:100,100),b=1.0*1:10000), by[:a], s=mminimum(_b, 1, window=10))) --> (10000,)
    @fact size(@select(darr(a=repmat(1:100,100),b=1.0*1:10000), by[:a], s=mmaximum(_b))) --> (100,)
    @fact size(@update(darr(a=repmat(1:100,100),b=1.0*1:10000), by[:a], s=mmaximum(_b))) --> (10000,)
    @fact size(@select(darr(a=repmat(1:100,100),b=1.0*1:10000), by[:a], s=mmedian(_b))) --> (100,)
    @fact size(@update(darr(a=repmat(1:100,100),b=1.0*1:10000), by[:a], s=mmedian(_b))) --> (10000,)
    @fact size(@select(darr(a=repmat(1:100,100),b=1.0*1:10000), by[:a], s=mmiddle(_b, 1, rev=true))) --> (100,)
    @fact size(@update(darr(a=repmat(1:100,100),b=1.0*1:10000), by[:a], s=mmiddle(_b))) --> (10000,)
    @fact size(@select(darr(a=repmat(1:100,100),b=1.0*1:10000), by[:a], s=mquantile(_b, 0.2, window=20))) --> (100,)
    @fact size(@update(darr(a=repmat(1:100,100),b=1.0*1:10000), by[:a], s=mquantile(_b, 0.5))) --> (10000,)
    @fact size(@select(darr(a=repmat(1:100,100),b=1.0*1:10000), by[:a], s=mquantile(_b, 0.1, rev=true))) --> (100,)
    @fact size(@update(darr(a=repmat(1:100,100),b=1.0*1:10000), by[:a], s=mquantile(_b, 0.8, rev=true))) --> (10000,)
  end
  context("misc tests") do
    @fact @select(larr(a=[1,1,2,3,4],b=enumeration([:a,:a,:b,:b,:b])), ct=length(_), by[:b]) --> larr(axis1=darr(b=[:a,:b]), ct=[2,3])
    @fact typeof(pickaxis(@select(larr(a=[1,1,2,3,4],b=enumeration([:a,:a,:b,:b,:b])), ct=length(_), by[:b]),1)) --> DataCubes.DictArray{Symbol,1,DataCubes.AbstractArrayWrapper{Nullable{Symbol},1,DataCubes.EnumerationArray{Symbol,1,Array{Int64,1},Int64}},Nullable{Symbol}}
    @fact @select(larr(a=[1,1,2,3,4],b=enumeration([1.0,1.0,2.0,1.0,2.0])), ct=length(_), by[:b]) --> @larr(axis1[b=[1.0,2.0]], ct=[3,2])
    @fact (@rap typeof pickaxis(_,1) @select(larr(a=[1,1,2,3,4],b=enumeration([1.0,1.0,2.0,1.0,2.0])), ct=length(_), by[:b])) --> DataCubes.DictArray{Symbol,1,DataCubes.AbstractArrayWrapper{Nullable{Float64},1,DataCubes.FloatNAArray{Float64,1,Array{Float64,1}}},Nullable{Float64}}
    @fact @select(larr(a=[1 2 3;4 5 6], axis2=darr(r=enumeration([:a,:a,:b]))), by[:r], s=sum(_a)) --> @larr(axis1[r=[:a,:b]], s=[12,9])
    @fact typeof(pickaxis(@select(larr(a=[1 2 3;4 5 6], axis2=darr(r=enumeration([:a,:a,:b]))), by[:r], s=sum(_a)),1)) --> DataCubes.DictArray{Symbol,1,DataCubes.AbstractArrayWrapper{Nullable{Symbol},1,DataCubes.EnumerationArray{Symbol,1,Array{Int64,1},Int64}},Nullable{Symbol}}
    @fact @select(larr(a=[1 2 3;4 5 6], axis2=darr(r=enumeration([:a,:a,:b]))), by[r=map(x->string(x.value),_r)], s=sum(_a)) --> @larr(axis1[r=["a","b"]], s=[12,9])
    @fact typeof(pickaxis(@select(larr(a=[1 2 3;4 5 6], axis2=darr(r=enumeration([:a,:a,:b]))), by[r=map(x->string(x.value),_r)], s=sum(_a)),1)) --> DataCubes.DictArray{Symbol,1,DataCubes.AbstractArrayWrapper{Nullable{String},1,Array{Nullable{String},1}},Nullable{String}}
    @fact @select(larr(a=[1 2 3;4 5 6], axis2=darr(r=[1.0,1.0,2.0])), by[:r], s=sum(_a)) --> @larr(axis1[r=[1.0,2.0]], s=[12,9])
    @fact typeof(pickaxis(@select(larr(a=[1 2 3;4 5 6], axis2=darr(r=[1.0,1.0,2.0])), by[:r], s=sum(_a)),1)) --> DataCubes.DictArray{Symbol,1,DataCubes.AbstractArrayWrapper{Nullable{Float64},1,DataCubes.FloatNAArray{Float64,1,Array{Float64,1}}},Nullable{Float64}}
    @fact @select(larr(a=[1 2 3;4 5 6], axis2=darr(r=[1.0,1.0,2.0])), by[r=broadcast(*, _r, 2)], s=sum(_a)) --> @larr(axis1[r=[2.0,4.0]], s=[12,9])
    @fact typeof(pickaxis(@select(larr(a=[1 2 3;4 5 6], axis2=darr(r=[1.0,1.0,2.0])), by[r=broadcast(*, _r, 2)], s=sum(_a)),1)) --> DataCubes.DictArray{Symbol,1,DataCubes.AbstractArrayWrapper{Nullable{Float64},1,DataCubes.FloatNAArray{Float64,1,Array{Float64,1}}},Nullable{Float64}}
    @fact typeof(view(@select(larr(a=rand(5,3),axis2=darr(r=enumeration([:X,:Y,:Z]))),:r),2,1:2).data.data[:r]) -->  typeof(view(@select(larr(a=rand(5,3),axis2=darr(r=enumeration([:X,:Y,:Z]))),:r),2,1:2).data.data[:r])
    @fact @select(@larr(a=[1 2 3;4 5 6], axis1[k=[:a,:b]], axis2[r=['X','Y','Z']]), by[:r], s=sum(_a)) --> larr(axis1=darr(r=['X','Y','Z']), s=[5,7,9])
    @fact @select(@larr(a=[1 2 3;4 5 6], axis1[k=[:a,:b]], axis2[r=['X','Y','Z']]), by[:r,:r1=>_r], s=sum(_a)) --> larr(axis1=darr(r=['X','Y','Z'],r1=['X','Y','Z']), s=[5,7,9])
    @fact @select(@larr(a=[1 2 3;4 5 6], axis1[k=[:a,:b]], axis2[r=['X','Y','Z']]), by[:r,:r1=>_r,r2=_r], s=sum(_a)) --> larr(axis1=darr(r=['X','Y','Z'],r1=['X','Y','Z'],r2=['X','Y','Z']), s=[5,7,9])
    @fact @select(@larr(a=[1 2 3;4 5 6], axis1[k=[:a,:b]], axis2[r=['X','Y','Z']]), by[:r,:r1=>_r,r2=_r,:r3=>_r], s=sum(_a)) --> larr(axis1=darr(r=['X','Y','Z'],r1=['X','Y','Z'],r2=['X','Y','Z'],r3=['X','Y','Z']), s=[5,7,9])
    @fact @select(@larr(a=[1 2 3;4 5 6], axis1[k=[:a,:b]], axis2[r=['X','Y','Z']]), by[:r,:r1=>_r,r2=_r,:r3=>_r,:r4=>_r], s=sum(_a)) --> larr(axis1=darr(r=['X','Y','Z'],r1=['X','Y','Z'],r2=['X','Y','Z'],r3=['X','Y','Z'],r4=['X','Y','Z']), s=[5,7,9])
    @fact @select(@larr(a=[1 2 3;4 5 6], axis1[k=[:a,:b]], axis2[r=['X','Y','Z']]), by[:r,:r1=>_r,r2=_r,:r3=>_r,:r4=>_r,r5=_r], s=sum(_a)) --> larr(axis1=darr(r=['X','Y','Z'],r1=['X','Y','Z'],r2=['X','Y','Z'],r3=['X','Y','Z'],r4=['X','Y','Z'],r5=['X','Y','Z']), s=[5,7,9])
    @fact @select(@larr(a=[1 2 3;4 5 6], axis1[k=[:a,:b]], axis2[r=['X','Y','Z']]), by[:r,:r1=>_r,r2=_r,:r3=>_r,:r4=>_r,r5=_r,r6=_r], s=sum(_a)) --> larr(axis1=darr(r=['X','Y','Z'],r1=['X','Y','Z'],r2=['X','Y','Z'],r3=['X','Y','Z'],r4=['X','Y','Z'],r5=['X','Y','Z'],r6=['X','Y','Z']), s=[5,7,9])
    @fact @select(@larr(a=[1 2 3;4 5 6], axis1[k=[:a,:b]], axis2[r=['X','Y','Z']]), by[:r,:r1=>_r,r2=_r,:r3=>_r,:r4=>_r,r5=_r,r6=_r,r7=_r], s=sum(_a)) --> larr(axis1=darr(r=['X','Y','Z'],r1=['X','Y','Z'],r2=['X','Y','Z'],r3=['X','Y','Z'],r4=['X','Y','Z'],r5=['X','Y','Z'],r6=['X','Y','Z'],r7=['X','Y','Z']), s=[5,7,9])
    @fact @update(darr(a=["a" "b" "cd";"efg" "h" "ij"]), b=length(_a)) --> darr(a=["a" "b" "cd";"efg" "h" "ij"], b=fill(6,2,3))
    @fact @update(larr(a=["a" "b" "cd";"efg" "h" "ij"],axis1=[:x,:y]), b=length(_a)) --> larr(a=["a" "b" "cd";"efg" "h" "ij"], b=fill(6,2,3), axis=[:x,:y])
    @fact size(@select(darr(a=repmat(1:100,100),b=1.0*1:10_000), by[:a], s=msum(_b))) --> (100,)
    @fact size(@update(darr(a=repmat(1:100,100),b=1.0*1:10_000), by[:a], s=msum(_b))) --> (10000,)
    @fact size(@select(darr(a=repmat(1:100,100),b=1.0*1:10000), by[:a], s=mprod(_b))) --> (100,)
    @fact size(@update(darr(a=repmat(1:100,100),b=1.0*1:10000), by[:a], s=mprod(_b, rev=true, window=10))) --> (10000,)
    @fact size(@select(darr(a=repmat(1:100,100),b=1.0*1:10000), by[:a], s=mmean(_b, window=10))) --> (100,)
    @fact size(@update(darr(a=repmat(1:100,100),b=1.0*1:10000), by[:a], s=mmean(_b, window=1))) --> (10000,)
    @fact size(@select(darr(a=repmat(1:100,100),b=1.0*1:10000), by[:a], s=mminimum(_b))) --> (100,)
    @fact size(@update(darr(a=repmat(1:100,100),b=1.0*1:10000), by[:a], s=mminimum(_b, 1, window=10))) --> (10000,)
    @fact size(@select(darr(a=repmat(1:100,100),b=1.0*1:10000), by[:a], s=mmaximum(_b))) --> (100,)
    @fact size(@update(darr(a=repmat(1:100,100),b=1.0*1:10000), by[:a], s=mmaximum(_b))) --> (10000,)
    @fact size(@select(darr(a=repmat(1:100,100),b=1.0*1:10000), by[:a], s=mmedian(_b))) --> (100,)
    @fact size(@update(darr(a=repmat(1:100,100),b=1.0*1:10000), by[:a], s=mmedian(_b))) --> (10000,)
    @fact size(@select(darr(a=repmat(1:100,100),b=1.0*1:10000), by[:a], s=mmiddle(_b, 1, rev=true))) --> (100,)
    @fact size(@update(darr(a=repmat(1:100,100),b=1.0*1:10000), by[:a], s=mmiddle(_b))) --> (10000,)
    @fact size(@select(darr(a=repmat(1:100,100),b=1.0*1:10000), by[:a], s=mquantile(_b, 0.2, window=20))) --> (100,)
    @fact size(@update(darr(a=repmat(1:100,100),b=1.0*1:10000), by[:a], s=mquantile(_b, 0.5))) --> (10000,)
    @fact size(@select(darr(a=repmat(1:100,100),b=1.0*1:10000), by[:a], s=mquantile(_b, 0.1, rev=true))) --> (100,)
    @fact size(@update(darr(a=repmat(1:100,100),b=1.0*1:10000), by[:a], s=mquantile(_b, 0.8, rev=true))) --> (10000,)
  end
end

end


text:
using TerminalClock
using Documenter

DocMeta.setdocmeta!(
    TerminalClock,
    :DocTestSetup,
    :(using TerminalClock);
    recursive=true,
)

makedocs(;
    modules=[TerminalClock],
    authors="Satoshi Terasaki <terasakisatoshi.math@gmail.com> and contributors",
    repo="https://github.com/AtelierArith/TerminalClock.jl/blob/{commit}{path}#{line}",
    sitename="TerminalClock.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://atelierarith.github.io/TerminalClock.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/AtelierArith/TerminalClock.jl",
    devbranch="main",
)

ext: jl


text:
using Logging

@testset "Test simple amplitude sampling functions" begin
    qubits = 5
    amplitudes = amplitudes_all(qubits)
    @test length(amplitudes) == 2^qubits
    @test all(length.(amplitudes) .== qubits)

    qubits = 5
    N = 10
    amplitudes = amplitudes_uniform(qubits, nothing, N)
    @test length(amplitudes) == N
end

@testset "Test simulation of small circuits" begin
    n = 5
    circ = create_ghz_circuit(n)

    results = with_logger(NullLogger()) do
        run_simulation(circ)
    end

    @test length(results) == 2^n
    @test results["0"^n] ≈ results["1"^n] ≈ 1/sqrt(2)
end


text:
using DiffEqBase, DiffEqJump
using Plots; plotlyjs()

tf      = .01
u0      = [200, 100, 150]

# MODEL SETUP

# DiffEqBiological Version:
using DiffEqBiological
rs = @reaction_network dtype begin
    k1, 2A --> B
    k2, B --> 2A
    k3, A + B --> C
    k4, C --> A + B
    k5, 3C --> 3A
end k1 k2 k3 k4 k5
rates = [1., 2., .5, .75, .25]

# model using mass action jumps
# ids: A = 1, B = 2, C = 3

# reactant stoichiometry is stored as a Vector{Vector{Pairs}}
# the kth entry in reactstoch is the stoichiometry of the kth reaction
reactstoch =
[
    [1 => 2],            # 2A
    [2 => 1],            # B
    [1 => 1, 2 => 1],    # A, B
    [3 => 1],            # C
    [3 => 3]             # 3C
]

# net stoichiometry is stored as a Vector{Vector{Pairs}}
# the kth entry in netstoch is the net stoichiometry of the kth reaction
netstoch =
[
    [1 => -2, 2 => 1],          # A -= 2, B += 1
    [1 => 2, 2 => -1],          # A += 2, B -= 1
    [1 => -1, 2 => -1, 3 => 1], # A -= 1, B -= 1, C += 1
    [1 => 1, 2 => 1, 3 => -1],  # A += 1, B += 1, C -= 1
    [1 => 3, 3 => -3]           # A += 3, C -= 3
]

# how to create the mass action jump type
majumps = MassActionJump(rates, reactstoch, netstoch)

prob = DiscreteProblem(u0, (0.0, tf), rates)

# solving with DiffEqBiological
jump_prob = JumpProblem(prob, Direct(), rs)
sol = solve(jump_prob, SSAStepper())
plothand = plot(sol, seriestype=:steppost, reuse=false)
display(plothand)

# solving with majumps
jump_prob = JumpProblem(prob, Direct(), majumps)
sol = solve(jump_prob, SSAStepper())
plothand = plot(sol, seriestype=:steppost, reuse=false)
display(plothand)

max_stars_repo_name: isaacsas/DiffEqJumpExtensions

max_issues_repo_name: isaacsas/DiffEqJumpExtensions

max_forks_repo_name: isaacsas/DiffEqJumpExtensions


text:
#####################################################
# Surya
# NREL
# November 2021
# Main function for EGRET JSON --> CSV files formatted according to SIIP Tabular Data Parser Requirements --> SIIP PSY System
#####################################################################################
function EGRET_TO_PSY(EGRET_json::Dict{String, Any};EGRET_json_RT::Union{Nothing, Dict{String, Any}} = nothing,location::Union{Nothing, String} = nothing,
                      ts_pointers_file::Union{Nothing, String} = nothing)

    if (ts_pointers_file === nothing)
    @warn "Time series pointers file type wasn't passed. Using timeseries_pointers.csv"
    ts_pointers_file = "CSV"
    end

    location, base_MVA,rt_flag =
    if (EGRET_json_RT !== nothing)
        parse_EGRET_JSON(EGRET_json,EGRET_json_RT=EGRET_json_RT)
    else
        parse_EGRET_JSON(EGRET_json)
    end

    if (rt_flag)
        sys_DA,sys_RT = parse_tabular_data(location,base_MVA,rt_flag,ts_pointers_file=ts_pointers_file)
        return sys_DA,sys_RT
    else
        sys_DA = parse_tabular_data(location,base_MVA,rt_flag,ts_pointers_file=ts_pointers_file)
        return sys_DA
    end
end


text:
module NumericalSoundSynthesis

include("musical-tunings.jl")
include("instruments.jl")
include("keyboard-piano.jl")

include("sine-instrument.jl")

end # module


text:
@testset "1.2.1 Quadratic" begin
    include("1.2.1.1 (a+b x+c x^2)^p.jl")
    include("1.2.1.2 (d+e x)^m (a+b x+c x^2)^p.jl")
    include("1.2.1.3 (d+e x)^m (f+g x) (a+b x+c x^2)^p.jl")
    include("1.2.1.4 (d+e x)^m (f+g x)^n (a+b x+c x^2)^p.jl")
    include("1.2.1.5 (a+b x+c x^2)^p (d+e x+f x^2)^q.jl")
    include("1.2.1.6 (g+h x)^m (a+b x+c x^2)^p (d+e x+f x^2)^q.jl")
    include("1.2.1.9 P(x) (d+e x)^m (a+b x+c x^2)^p.jl")
end

max_stars_repo_path: test/1 Algebraic functions/1.2 Trinomial products/1.2.1 Quadratic/.jl


text:
"""
    GeneralizedLinearMixedModel

Generalized linear mixed-effects model representation

# Fields
- `LMM`: a [`LinearMixedModel`](@ref) - the local approximation to the GLMM.
- `β`: the pivoted and possibly truncated fixed-effects vector
- `β₀`: similar to `β`. Used in the PIRLS algorithm if step-halving is needed.
- `θ`: covariance parameter vector
- `b`: similar to `u`, equivalent to `broadcast!(*, b, LMM.Λ, u)`
- `u`: a vector of matrices of random effects
- `u₀`: similar to `u`.  Used in the PIRLS algorithm if step-halving is needed.
- `resp`: a `GlmResp` object
- `η`: the linear predictor
- `wt`: vector of prior case weights, a value of `T[]` indicates equal weights.
The following fields are used in adaptive Gauss-Hermite quadrature, which applies
only to models with a single random-effects term, in which case their lengths are
the number of levels in the grouping factor for that term.  Otherwise they are
zero-length vectors.
- `devc`: vector of deviance components
- `devc0`: vector of deviance components at offset of zero
- `sd`: approximate standard deviation of the conditional density
- `mult`: multiplier

# Properties

In addition to the fieldnames, the following names are also accessible through the `.` extractor

- `theta`: synonym for `θ`
- `beta`: synonym for `β`
- `σ` or `sigma`: common scale parameter (value is `NaN` for distributions without a scale parameter)
- `lowerbd`: vector of lower bounds on the combined elements of `β` and `θ`
- `formula`, `trms`, `A`, `L`, and `optsum`: fields of the `LMM` field
- `X`: fixed-effects model matrix
- `y`: response vector

"""
struct GeneralizedLinearMixedModel{T<:AbstractFloat,D<:Distribution} <: MixedModel{T}
    LMM::LinearMixedModel{T}
    β::Vector{T}
    β₀::Vector{T}
    θ::Vector{T}
    b::Vector{Matrix{T}}
    u::Vector{Matrix{T}}
    u₀::Vector{Matrix{T}}
    resp::GLM.GlmResp
    η::Vector{T}
    wt::Vector{T}
    devc::Vector{T}
    devc0::Vector{T}
    sd::Vector{T}
    mult::Vector{T}
end

function StatsBase.coef(m::GeneralizedLinearMixedModel{T}) where {T}
    piv = m.LMM.feterm.piv
    return invpermute!(copyto!(fill(T(-0.0), length(piv)), m.β), piv)
end

function StatsBase.coeftable(m::GeneralizedLinearMixedModel)
    co = coef(m)
    se = stderror(m)
    z = co ./ se
    pvalue = ccdf.(Chisq(1), abs2.(z))
    return CoefTable(
        hcat(co, se, z, pvalue),
        ["Coef.", "Std. Error", "z", "Pr(>|z|)"],
        coefnames(m),
        4, # pvalcol
        3, # teststatcol
    )
end

"""
    deviance(m::GeneralizedLinearMixedModel{T}, nAGQ=1)::T where {T}

Return the deviance of `m` evaluated by the Laplace approximation (`nAGQ=1`)
or `nAGQ`-point adaptive Gauss-Hermite quadrature.

If the distribution `D` does not have a scale parameter the Laplace approximation
is the squared length of the conditional modes, ``u``, plus the determinant
of ``Λ'Z'WZΛ + I``, plus the sum of the squared deviance residuals.
"""
function StatsBase.deviance(m::GeneralizedLinearMixedModel{T}, nAGQ=1) where {T}
    nAGQ == 1 && return T(sum(m.resp.devresid) + logdet(m) + sum(u -> sum(abs2, u), m.u))
    u = vec(first(m.u))
    u₀ = vec(first(m.u₀))
    copyto!(u₀, u)
    ra = RaggedArray(m.resp.devresid, first(m.LMM.reterms).refs)
    devc0 = sum!(map!(abs2, m.devc0, u), ra)  # the deviance components at z = 0
    sd = map!(inv, m.sd, first(m.LMM.L).diag)
    mult = fill!(m.mult, 0)
    devc = m.devc
    for (z, w) in GHnorm(nAGQ)
        if !iszero(w)
            if iszero(z)  # devc == devc0 in this case
                mult .+= w
            else
                @. u = u₀ + z * sd
                updateη!(m)
                sum!(map!(abs2, devc, u), ra)
                @. mult += exp((abs2(z) + devc0 - devc) / 2) * w
            end
        end
    end
    copyto!(u, u₀)
    updateη!(m)
    return sum(devc0) - 2 * (sum(log, mult) + sum(log, sd))
end

StatsBase.deviance(m::GeneralizedLinearMixedModel) = deviance(m, m.optsum.nAGQ)

fixef(m::GeneralizedLinearMixedModel) = m.β

function fixef!(v::AbstractVector{T}, m::GeneralizedLinearMixedModel{T}) where {T}
    return copyto!(fill!(v, -zero(T)), m.β)
end

objective(m::GeneralizedLinearMixedModel) = deviance(m)

"""
    GLM.wrkresp!(v::AbstractVector{T}, resp::GLM.GlmResp{AbstractVector{T}})

A copy of a method from GLM that generalizes the types in the signature
"""
function GLM.wrkresp!(
    v::AbstractVector{T}, r::GLM.GlmResp{Vector{T}}
) where {T<:AbstractFloat}
    v .= r.eta .+ r.wrkresid
    isempty(r.offset) && return v
    return v .-= r.offset
end

"""
    deviance!(m::GeneralizedLinearMixedModel, nAGQ=1)

Update `m.η`, `m.μ`, etc., install the working response and working weights in
`m.LMM`, update `m.LMM.A` and `m.LMM.R`, then evaluate the `deviance`.
"""
function deviance!(m::GeneralizedLinearMixedModel, nAGQ=1)
    updateη!(m)
    GLM.wrkresp!(m.LMM.y, m.resp)
    reweight!(m.LMM, m.resp.wrkwt)
    return deviance(m, nAGQ)
end

function GLM.dispersion(m::GeneralizedLinearMixedModel{T}, sqr::Bool=false) where {T}
    # adapted from GLM.dispersion(::AbstractGLM, ::Bool)
    # TODO: PR for a GLM.dispersion(resp::GLM.GlmResp, dof_residual::Int, sqr::Bool)
    r = m.resp
    if dispersion_parameter(r.d)
        s = sum(wt * abs2(re) for (wt, re) in zip(r.wrkwt, r.wrkresid)) / dof_residual(m)
        sqr ? s : sqrt(s)
    else
        one(T)
    end
end

GLM.dispersion_parameter(m::GeneralizedLinearMixedModel) = dispersion_parameter(m.resp.d)

Distributions.Distribution(m::GeneralizedLinearMixedModel{T,D}) where {T,D} = D

function fit(
    ::Type{GeneralizedLinearMixedModel},
    f::FormulaTerm,
    tbl,
    d::Distribution=Normal(),
    l::Link=canonicallink(d);
    wts=[],
    contrasts=Dict{Symbol,Any}(),
    offset=[],
    verbose::Bool=false,
    fast::Bool=false,
    nAGQ::Integer=1,
    progress::Bool=true,
    thin::Int=typemax(Int),
)
    return fit(
        GeneralizedLinearMixedModel,
        f,
        columntable(tbl),
        d,
        l;
        wts,
        offset,
        contrasts,
        verbose,
        fast,
        nAGQ,
        progress,
        thin,
    )
end

function fit(
    ::Type{GeneralizedLinearMixedModel},
    f::FormulaTerm,
    tbl::Tables.ColumnTable,
    d::Distribution,
    l::Link=canonicallink(d);
    wts=[],
    contrasts=Dict{Symbol,Any}(),
    offset=[],
    verbose::Bool=false,
    fast::Bool=false,
    nAGQ::Integer=1,
    progress::Bool=true,
    thin::Int=typemax(Int),
)
    return fit!(
        GeneralizedLinearMixedModel(
            f, tbl, d, l; wts=wts, offset=offset, contrasts=contrasts
        );
        verbose,
        fast,
        nAGQ,
        progress,
        thin,
    )
end

function fit(
    ::Type{MixedModel},
    f::FormulaTerm,
    tbl,
    d::Distribution,
    l::Link=canonicallink(d);
    wts=[],
    contrasts=Dict{Symbol,Any}(),
    offset=[],
    verbose::Bool=false,
    REML::Bool=false,
    fast::Bool=false,
    nAGQ::Integer=1,
    progress::Bool=true,
    thin::Int=typemax(Int),
)
    return fit(
        GeneralizedLinearMixedModel,
        f,
        tbl,
        d,
        l;
        wts,
        contrasts,
        offset,
        verbose,
        fast,
        nAGQ,
        progress,
        thin,
    )
end

"""
    fit!(m::GeneralizedLinearMixedModel; fast=false, nAGQ=1,
                                         verbose=false, progress=true,
                                         thin::Int=1)

Optimize the objective function for `m`.

When `fast` is `true` a potentially much faster but slightly less accurate algorithm, in
which `pirls!` optimizes both the random effects and the fixed-effects parameters,
is used.

If `progress` is `true`, the default, a `ProgressMeter.ProgressUnknown` counter is displayed.
during the iterations to minimize the deviance.  There is a delay before this display is initialized
and it may not be shown at all for models that are optimized quickly.

If `verbose` is `true`, then both the intermediate results of both the nonlinear optimization and PIRLS are also displayed on standard output.

At every `thin`th iteration  is recorded in `fitlog`, optimization progress is saved in `m.optsum.fitlog`.
"""
function fit!(
    m::GeneralizedLinearMixedModel{T};
    verbose::Bool=false,
    fast::Bool=false,
    nAGQ::Integer=1,
    progress::Bool=true,
    thin::Int=typemax(Int),
) where {T}
    β = m.β
    lm = m.LMM
    optsum = lm.optsum

    if optsum.feval > 0
        throw(ArgumentError("This model has already been fitted. Use refit!() instead."))
    end

    if !fast
        optsum.lowerbd = vcat(fill!(similar(β), T(-Inf)), optsum.lowerbd)
        optsum.initial = vcat(β, m.θ)
        optsum.final = copy(optsum.initial)
    end
    setpar! = fast ? setθ! : setβθ!
    prog = ProgressUnknown("Minimizing"; showspeed=true)
    # start from zero for the initial call to obj before optimization
    iter = 0
    fitlog = optsum.fitlog
    function obj(x, g)
        isempty(g) || throw(ArgumentError("g should be empty for this objective"))
        val = deviance(pirls!(setpar!(m, x), fast, verbose), nAGQ)
        iszero(rem(iter, thin)) && push!(fitlog, (copy(x), val))
        verbose && println(round(val; digits=5), " ", x)
        progress && ProgressMeter.next!(prog; showvalues=[(:objective, val)])
        iter += 1
        return val
    end
    opt = Opt(optsum)
    NLopt.min_objective!(opt, obj)
    optsum.finitial = obj(optsum.initial, T[])
    empty!(fitlog)
    push!(fitlog, (copy(optsum.initial), optsum.finitial))
    fmin, xmin, ret = NLopt.optimize(opt, copyto!(optsum.final, optsum.initial))
    ProgressMeter.finish!(prog)
    ## check if very small parameter values bounded below by zero can be set to zero
    xmin_ = copy(xmin)
    for i in eachindex(xmin_)
        if iszero(optsum.lowerbd[i]) && zero(T) < xmin_[i] < T(0.001)
            xmin_[i] = zero(T)
        end
    end
    loglength = length(fitlog)
    if xmin ≠ xmin_
        if (zeroobj = obj(xmin_, T[])) ≤ (fmin + 1.e-5)
            fmin = zeroobj
            copyto!(xmin, xmin_)
        elseif length(fitlog) > loglength
            # remove unused extra log entry
            pop!(fitlog)
        end
    end
    ## ensure that the parameter values saved in m are xmin
    pirls!(setpar!(m, xmin), fast, verbose)
    optsum.nAGQ = nAGQ
    optsum.feval = opt.numevals
    optsum.final = xmin
    optsum.fmin = fmin
    optsum.returnvalue = ret
    _check_nlopt_return(ret)
    return m
end

StatsBase.fitted(m::GeneralizedLinearMixedModel) = m.resp.mu

function GeneralizedLinearMixedModel(
    f::FormulaTerm,
    tbl,
    d::Distribution,
    l::Link=canonicallink(d);
    wts=[],
    offset=[],
    contrasts=Dict{Symbol,Any}(),
)
    return GeneralizedLinearMixedModel(
        f, Tables.columntable(tbl), d, l; wts=wts, offset=offset, contrasts=contrasts
    )
end

function GeneralizedLinearMixedModel(
    f::FormulaTerm,
    tbl::Tables.ColumnTable,
    d::Normal,
    l::IdentityLink;
    wts=[],
    offset=[],
    contrasts=Dict{Symbol,Any}(),
)
    return throw(
        ArgumentError("use LinearMixedModel for Normal distribution with IdentityLink")
    )
end

function GeneralizedLinearMixedModel(
    f::FormulaTerm,
    tbl::Tables.ColumnTable,
    d::Distribution,
    l::Link=canonicallink(d);
    wts=[],
    offset=[],
    contrasts=Dict{Symbol,Any}(),
)
    if isa(d, Binomial) && isempty(wts)
        d = Bernoulli()
    end
    (isa(d, Normal) && isa(l, IdentityLink)) && throw(
        ArgumentError("use LinearMixedModel for Normal distribution with IdentityLink")
    )

    if !any(isa(d, dist) for dist in (Bernoulli, Binomial, Poisson))
        @warn """Results for families with a dispersion parameter are not reliable.
                 It is best to avoid trying to fit such models in MixedModels until
                 the authors gain a better understanding of those cases."""
    end

    LMM = LinearMixedModel(f, tbl; contrasts=contrasts, wts=wts)
    y = copy(LMM.y)
    # the sqrtwts field must be the correct length and type but we don't know those
    # until after the model is constructed if wt is empty.  Because a LinearMixedModel
    # type is immutable, another one must be created.
    if isempty(wts)
        LMM = LinearMixedModel(
            LMM.formula,
            LMM.reterms,
            LMM.Xymat,
            LMM.feterm,
            fill!(similar(y), 1),
            LMM.parmap,
            LMM.dims,
            LMM.A,
            LMM.L,
            LMM.optsum,
        )
    end
    updateL!(LMM)
    # fit a glm to the fixed-effects only
    T = eltype(LMM.Xymat)
    gl = glm(LMM.X, y, d, l; wts=convert(Vector{T}, wts), offset=convert(Vector{T}, offset))
    β = coef(gl)
    u = [fill(zero(eltype(y)), vsize(t), nlevs(t)) for t in LMM.reterms]
    # vv is a template vector used to initialize fields for AGQ
    # it is empty unless there is a single random-effects term
    vv = length(u) == 1 ? vec(first(u)) : similar(y, 0)

    res = GeneralizedLinearMixedModel{T,typeof(d)}(
        LMM,
        β,
        copy(β),
        LMM.θ,
        copy.(u),
        u,
        zero.(u),
        gl.rr,
        similar(y),
        oftype(y, wts),
        similar(vv),
        similar(vv),
        similar(vv),
        similar(vv),
    )
    deviance!(res, 1)
    return res
end

function Base.getproperty(m::GeneralizedLinearMixedModel, s::Symbol)
    if s == :theta
        m.θ
    elseif s == :coef
        coef(m)
    elseif s == :beta
        m.β
    elseif s == :objective
        objective(m)
    elseif s ∈ (:σ, :sigma)
        sdest(m)
    elseif s == :σs
        σs(m)
    elseif s == :σρs
        σρs(m)
    elseif s ∈ (:A, :L, :optsum, :reterms, :Xymat, :feterm, :formula, :parmap)
        getfield(m.LMM, s)
    elseif s ∈ (:dims, :λ, :lowerbd, :corr, :PCA, :rePCA, :X)
        getproperty(m.LMM, s)
    elseif s == :y
        m.resp.y
    else
        getfield(m, s)
    end
end

# this copy behavior matches the implicit copy behavior
# for LinearMixedModel. So this is then different than m.θ,
# which returns a reference to the same array
getθ(m::GeneralizedLinearMixedModel) = copy(m.θ)
getθ!(v::AbstractVector{T}, m::GeneralizedLinearMixedModel{T}) where {T} = copyto!(v, m.θ)

StatsBase.islinear(m::GeneralizedLinearMixedModel) = isa(GLM.Link, GLM.IdentityLink)

GLM.Link(m::GeneralizedLinearMixedModel) = GLM.Link(m.resp)

function StatsBase.loglikelihood(m::GeneralizedLinearMixedModel{T}) where {T}
    accum = zero(T)
    # adapted from GLM.jl
    # note the use of loglik_obs to handle the different parameterizations
    # of various response distributions which may not just be location+scale
    r = m.resp
    wts = r.wts
    y = r.y
    mu = r.mu
    d = r.d
    if length(wts) == length(y)
        ϕ = deviance(r) / sum(wts)
        @inbounds for i in eachindex(y, mu, wts)
            accum += GLM.loglik_obs(d, y[i], mu[i], wts[i], ϕ)
        end
    else
        ϕ = deviance(r) / length(y)
        @inbounds for i in eachindex(y, mu)
            accum += GLM.loglik_obs(d, y[i], mu[i], 1, ϕ)
        end
    end
    return accum - (mapreduce(u -> sum(abs2, u), +, m.u) + logdet(m)) / 2
end

function Base.propertynames(m::GeneralizedLinearMixedModel, private::Bool=false)
    return (
        :A,
        :L,
        :theta,
        :beta,
        :coef,
        :λ,
        :σ,
        :sigma,
        :X,
        :y,
        :lowerbd,
        :objective,
        :σρs,
        :σs,
        :corr,
        :vcov,
        :PCA,
        :rePCA,
        (
            if private
                fieldnames(GeneralizedLinearMixedModel)
            else
                (:LMM, :β, :θ, :b, :u, :resp, :wt)
            end
        )...,
    )
end

"""
    pirls!(m::GeneralizedLinearMixedModel)

Use Penalized Iteratively Reweighted Least Squares (PIRLS) to determine the conditional
modes of the random effects.

When `varyβ` is true both `u` and `β` are optimized with PIRLS.  Otherwise only `u` is
optimized and `β` is held fixed.

Passing `verbose = true` provides verbose output of the iterations.
"""
function pirls!(
    m::GeneralizedLinearMixedModel{T}, varyβ=false, verbose=false; maxiter::Integer=10
) where {T}
    u₀ = m.u₀
    u = m.u
    β = m.β
    β₀ = m.β₀
    lm = m.LMM
    for j in eachindex(u)         # start from u all zeros
        copyto!(u₀[j], fill!(u[j], 0))
    end
    if varyβ
        copyto!(β₀, β)
        Llast = last(lm.L)
        pp1 = size(Llast, 1)
        Ltru = view(Llast, pp1, 1:(pp1 - 1)) # name read as L'u
    end
    obj₀ = deviance!(m) * 1.0001
    if verbose
        print("varyβ = ", varyβ, ", obj₀ = ", obj₀)
        if varyβ
            print(", β = ")
            show(β)
        end
        println()
    end
    for iter in 1:maxiter
        varyβ && ldiv!(adjoint(feL(m)), copyto!(β, Ltru))
        ranef!(u, m.LMM, β, true) # solve for new values of u
        obj = deviance!(m)        # update GLM vecs and evaluate Laplace approx
        verbose && println(lpad(iter, 4), ": ", obj)
        nhalf = 0
        while obj > obj₀
            nhalf += 1
            if nhalf > 10
                if iter < 2
                    throw(ErrorException("number of averaging steps > 10"))
                end
                break
            end
            for i in eachindex(u)
                map!(average, u[i], u[i], u₀[i])
            end
            varyβ && map!(average, β, β, β₀)
            obj = deviance!(m)
            verbose && println(lpad(nhalf, 8), ", ", obj)
        end
        if isapprox(obj, obj₀; atol=0.00001)
            break
        end
        copyto!.(u₀, u)
        copyto!(β₀, β)
        obj₀ = obj
    end
    return m
end

ranef(m::GeneralizedLinearMixedModel; uscale::Bool=false) = ranef(m.LMM; uscale=uscale)

LinearAlgebra.rank(m::GeneralizedLinearMixedModel) = m.LMM.feterm.rank

"""
    refit!(m::GeneralizedLinearMixedModel[, y::Vector];
           fast::Bool = (length(m.θ) == length(m.optsum.final)),
           nAGQ::Integer = m.optsum.nAGQ,
           kwargs...)

Refit the model `m` after installing response `y`.

If `y` is omitted the current response vector is used.

If not specified, the `fast` and `nAGQ` options from the previous fit are used.
`kwargs` are the same as [`fit!`](@ref)
"""
function refit!(
    m::GeneralizedLinearMixedModel;
    fast::Bool=(length(m.θ) == length(m.optsum.final)),
    nAGQ::Integer=m.optsum.nAGQ,
    kwargs...,
)
    return fit!(unfit!(m); fast=fast, nAGQ=nAGQ, kwargs...)
end

function refit!(m::GeneralizedLinearMixedModel, y; kwargs...)
    m_resp_y = m.resp.y
    length(y) == size(m_resp_y, 1) || throw(DimensionMismatch(""))
    copyto!(m_resp_y, y)
    return refit!(m; kwargs...)
end

"""
    setβθ!(m::GeneralizedLinearMixedModel, v)

Set the parameter vector, `:βθ`, of `m` to `v`.

`βθ` is the concatenation of the fixed-effects, `β`, and the covariance parameter, `θ`.
"""
function setβθ!(m::GeneralizedLinearMixedModel, v)
    setβ!(m, v)
    return setθ!(m, view(v, (length(m.β) + 1):length(v)))
end

function setβ!(m::GeneralizedLinearMixedModel, v)
    β = m.β
    copyto!(β, view(v, 1:length(β)))
    return m
end

function setθ!(m::GeneralizedLinearMixedModel, v)
    setθ!(m.LMM, copyto!(m.θ, v))
    return m
end

function Base.setproperty!(m::GeneralizedLinearMixedModel, s::Symbol, y)
    if s == :β
        setβ!(m, y)
    elseif s == :θ
        setθ!(m, y)
    elseif s == :βθ
        setβθ!(m, y)
    else
        setfield!(m, s, y)
    end
end

"""
    sdest(m::GeneralizedLinearMixedModel)

Return the estimate of the dispersion, i.e. the standard deviation of the per-observation noise.

For models with a dispersion parameter ϕ, this is simply ϕ. For models without a
dispersion parameter, this value is `missing`. This differs from `disperion`,
which returns `1` for models without a dispersion parameter.

For Gaussian models, this parameter is often called σ.
"""
function sdest(m::GeneralizedLinearMixedModel{T}) where {T}
    return dispersion_parameter(m) ? dispersion(m, false) : missing
end

function Base.show(
    io::IO, ::MIME"text/plain", m::GeneralizedLinearMixedModel{T,D}
) where {T,D}
    if m.optsum.feval < 0
        @warn("Model has not been fit")
        return nothing
    end
    nAGQ = m.LMM.optsum.nAGQ
    println(io, "Generalized Linear Mixed Model fit by maximum likelihood (nAGQ = $nAGQ)")
    println(io, "  ", m.LMM.formula)
    println(io, "  Distribution: ", D)
    println(io, "  Link: ", Link(m), "\n")
    nums = Ryu.writefixed.([loglikelihood(m), deviance(m), aic(m), aicc(m), bic(m)], 4)
    fieldwd = max(maximum(textwidth.(nums)) + 1, 11)
    for label in [" logLik", " deviance", "AIC", "AICc", "BIC"]
        print(io, rpad(lpad(label, (fieldwd + textwidth(label)) >> 1), fieldwd))
    end
    println(io)
    print.(Ref(io), lpad.(nums, fieldwd))
    println(io)
    println(io)

    show(io, VarCorr(m))

    print(io, " Number of obs: $(length(m.y)); levels of grouping factors: ")
    join(io, nlevs.(m.reterms), ", ")
    println(io)

    println(io, "\nFixed-effects parameters:")
    return show(io, coeftable(m))
end

Base.show(io::IO, m::GeneralizedLinearMixedModel) = show(io, MIME"text/plain"(), m)

function stderror!(v::AbstractVector{T}, m::GeneralizedLinearMixedModel{T}) where {T}
    # initialize to appropriate NaN for rank-deficient case
    fill!(v, zero(T) / zero(T))

    # the inverse permutation is done here.
    # if this is changed to access the permuted
    # model matrix directly, then don't forget to add
    # in the inverse permutation
    vcovmat = vcov(m)

    for idx in 1:size(vcovmat, 1)
        v[idx] = sqrt(vcovmat[idx, idx])
    end

    return v
end

function unfit!(model::GeneralizedLinearMixedModel{T}) where {T}
    deviance!(model, 1)
    reevaluateAend!(model.LMM)

    reterms = model.LMM.reterms
    optsum = model.LMM.optsum
    # we need to reset optsum so that it
    # plays nice with the modifications fit!() does
    optsum.lowerbd = mapfoldl(lowerbd, vcat, reterms)
    optsum.initial = mapfoldl(getθ, vcat, reterms)
    optsum.final = copy(optsum.initial)
    optsum.xtol_abs = fill!(copy(optsum.initial), 1.0e-10)
    optsum.initial_step = T[]
    optsum.feval = -1

    return model
end

"""
    updateη!(m::GeneralizedLinearMixedModel)

Update the linear predictor, `m.η`, from the offset and the `B`-scale random effects.
"""
function updateη!(m::GeneralizedLinearMixedModel{T}) where {T}
    η = m.η
    b = m.b
    u = m.u
    reterms = m.LMM.reterms
    mul!(η, modelmatrix(m), m.β)
    for i in eachindex(b)
        mul!(η, reterms[i], vec(mul!(b[i], reterms[i].λ, u[i])), one(T), one(T))
    end
    GLM.updateμ!(m.resp, η)
    return m
end

"""
    varest(m::GeneralizedLinearMixedModel)

Returns the estimate of ϕ², the variance of the conditional distribution of Y given B.

For models with a dispersion parameter ϕ, this is simply ϕ². For models without a
dispersion parameter, this value is `missing`. This differs from `disperion`,
which returns `1` for models without a dispersion parameter.

For Gaussian models, this parameter is often called σ².
"""
function varest(m::GeneralizedLinearMixedModel{T}) where {T}
    return dispersion_parameter(m) ? dispersion(m, true) : missing
end

function StatsBase.weights(m::GeneralizedLinearMixedModel{T}) where {T}
    wts = m.wt
    return isempty(wts) ? ones(T, nobs(m)) : wts
end

# delegate GLMM method to LMM field
for f in (:feL, :fetrm, :fixefnames, :(LinearAlgebra.logdet), :lowerbd, :PCA, :rePCA)
    @eval begin
        $f(m::GeneralizedLinearMixedModel) = $f(m.LMM)
    end
end


text:
# coding: utf-8

# In[1]:

#fractal is one of the interesting topics in geometry
#it is usually described by a recursive function
#voila,here we are!
using Plots

# In[2]:

#define starting and ending points
#and the order of hilbert curve
n=4
starting_point=(0,0)
ending_point=(2^n-1,0)

# In[3]:

#the code for hilbert curve may look intimidating at the first glance
#it is extremely simple and straight forward
#julia version is slightly different
#as we cannot use yield like python
#coroutine in julia is too complicated
#maybe another day...
function hilbert_curve(coordinates,point1,point2,n)

    #unpack
    x_start,y_start=point1
    x_end,y_end=point2

    #the function consists four different parts
    #which are 4 different rotations plus flips of a second order hilbert curve
    #0 degree second order hilbert curve
    if x_start<x_end && y_start==y_end

        #compute the grid length
        L=x_end-x_start

        #keypoints are starting and ending points
        #of 4 different first order hilbert curves
        #on a second order hilbert curve
        keypoints=[(x_start,y_start),(x_start,y_start-(L-1)/2),
        (x_start,y_start-(L-1)/2-1),(x_start+(L-1)/2,y_start-(L-1)/2-1),
        (x_start+(L-1)/2+1,y_start-(L-1)/2-1),(x_start+L,y_start-(L-1)/2-1),
        (x_start+L,y_start-(L-1)/2),(x_start+L,y_start)]

        #base case
        if n==1
            append!(coordinates,keypoints)
        else

            #recursion
            for i in 1:2:8
                hilbert_curve(coordinates,
                    keypoints[i],keypoints[i+1],n-1)
            end
        end
    end

    #180 degree second order hilbert curve
    if x_start>x_end && y_start==y_end
        L=x_start-x_end
        keypoints=[(x_start,y_start),(x_start,y_start+(L-1)/2),
        (x_start,y_start+(L-1)/2+1),(x_start-(L-1)/2,y_start+(L-1)/2+1),
        (x_start-(L-1)/2-1,y_start+(L-1)/2+1),(x_start-L,y_start+(L-1)/2+1),
        (x_start-L,y_start+(L-1)/2),(x_start-L,y_start)]
        if n==1
            append!(coordinates,keypoints)
        else
            for i in 1:2:8
                hilbert_curve(coordinates,
                    keypoints[i],keypoints[i+1],n-1)
            end
        end
    end

    #clockwise 90 degree horizontal flipped second order hilbert curve
    if x_start==x_end && y_start>y_end
        L=y_start-y_end
        keypoints=[(x_start,y_start),(x_start+(L-1)/2,y_start),
        (x_start+(L-1)/2+1,y_start),(x_start+(L-1)/2+1,y_start-(L-1)/2),
        (x_start+(L-1)/2+1,y_start-(L-1)/2-1),(x_start+(L-1)/2+1,y_start-L),
        (x_start+(L-1)/2,y_start-L),(x_start,y_start-L)]
        if n==1
            append!(coordinates,keypoints)
        else
            for i in 1:2:8
                hilbert_curve(coordinates,
                    keypoints[i],keypoints[i+1],n-1)
            end
        end
    end

    #clockwise 270 degree horizontal flipped second order hilbert curve
    if x_start==x_end && y_start<y_end
        L=y_end-y_start
        keypoints=[(x_start,y_start),(x_start-(L-1)/2,y_start),
        (x_start-(L-1)/2-1,y_start),(x_start-(L-1)/2-1,y_start+(L-1)/2),
        (x_start-(L-1)/2-1,y_start+(L-1)/2+1),(x_start-(L-1)/2-1,y_start+L),
        (x_start-(L-1)/2,y_start+L),(x_start,y_start+L)]
        if n==1
            append!(coordinates,keypoints)
        else
            for i in 1:2:8
                hilbert_curve(coordinates,
                    keypoints[i],keypoints[i+1],n-1)
            end
        end
    end
    coordinates
end

# In[4]:

#get coordinates
coordinates=hilbert_curve([],
    starting_point,ending_point,n);

# In[5]:

#viz
gr(size=(250,250))
fig=plot(legend=false,grid=false,axis=false,ticks=false)
plot!([i[1] for i in coordinates],
        [i[2] for i in coordinates])
fig

max_stars_repo_path: hilbert curve.jl

max_issues_repo_path: hilbert curve.jl

max_forks_repo_path: hilbert curve.jl


text:
using HypothesisTests
using Test
using StatsBase
using DelimitedFiles

@testset "Bartlett's test" begin
    # Columns are type, length, left, right, bottom, top, diag
    swiss = readdlm(joinpath(@__DIR__, "data", "swiss3.txt"))
    genuine = convert(Matrix{Float64}, swiss[view(swiss, :, 1) .== "real", 2:end])
    counterfeit = convert(Matrix{Float64}, swiss[view(swiss, :, 1) .== "fake", 2:end])

    b = BartlettTest(genuine, counterfeit)
    @test nobs(b) == (100, 100)
    @test dof(b) == 21
    @test pvalue(b) ≈ 0.0 atol=1e-10
    @test b.L′ ≈ 121.8991235 atol=1e-6
    @test occursin("reject h_0", sprint(show, b))
end


text:
using Printf

# Note that like enums, structs must be declared outside of functions
struct customer
    name::String
    balance::Float32
    age::Int
end

mutable struct mutable_customer
    name::String
    balance::Float32
    age::Int
end
function struct_examples()
    println("--------------------------------------")
    println("Struct examples")

    println()

    println("We can create instances of structs and set their values and then get them again")
    bob = customer("Bob Smith", 10.50, 46)
    println(bob)
    println(bob.name)

    #Note that structs are immutable and hence the following line would cause an error!
    #bob.name = "Robert"
    println("If we want to change values, we have to create a mutable struct")
    james = mutable_customer("James", 1234, 21)
    println(james)
    println(james.name)
    println("...now we can change the name")
    james.name = "Jim-bob"
    println(james)
    println(james.name)
end

max_stars_repo_name: James-P-D/JuliaDump

max_issues_repo_name: James-P-D/JuliaDump

max_forks_repo_name: James-P-D/JuliaDump


text:
using Documenter, Robust32s

makedocs(
    modules = [Robust32s],
    sitename = "Robust32s",
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    pages  = Any[
        "Overview" => "index.md",
        "How It Works" => "howitworks.md",
        "Helpful Hints" => "helpfulhints.md"
    ]
)

deploydocs(
    repo = "github.com/JeffreySarnoff/Robust32s.jl.git",
    target = "build"
)


text:
#Wigner 6j-symbol calculator.
#Translated and modified from function by Amita B Deb, Clarendon Lab. 2007.

# Calculates { j1, j2 ,j3}  using Racah formula. See: Sobelman: Atomic Spectra and Radiative Transitions.
#              J1  J2  J3

function Wigner6j(j1,j2,j3,J1,J2,J3)

  check = Wigner6jcheck(j1,j2,j3,J1,J2,J3)

  if check==1
    #Finding triangular coefficients
    tri1 = triangle_coeff(j1,j2,j3)
    tri2 = triangle_coeff(j1,J2,J3)
    tri3 = triangle_coeff(J1,j2,J3)
    tri4 = triangle_coeff(J1,J2,j3)

    #Finding the range of summation in Racah formula
    a = zeros(4)
    a[1] = j1 + j2 + j3
    a[2] = j1 + J2 + J3
    a[3] = J1 + j2 + J3
    a[4] = J1 + J2 + j3

    rangei = maximum(a)

    k = zeros(3)
    k[1] = j1 + j2 + J1 + J2
    k[2] = j2 + j3 + J2 + J3
    k[3] = j3 + j1 + J3 + J1

    rangef = minimum(k)

    Wignerhold = 0

    for tt = rangei:rangef
      Wignerhold = Wignerhold+ (tri1*tri2*tri3*tri4)^(0.5)*((-1)^tt)*factorial(tt+1)/fung(tt,j1,j2,j3,J1,J2,J3)
    end
    return Wigner = Wignerhold
  else
    return Wigner = 0
  end
end

max_stars_repo_name: CSChisholm/Rydberg

max_issues_repo_name: CSChisholm/Rydberg

max_forks_repo_name: CSChisholm/Rydberg


text:
using DiffEqBase, DiffEqJump
using Test

reactstoch =
[
    [1 => 1]
]

netstoch =
[
    [1 => -1]
]

rates = [1.]
u0 = [10]
spec_to_dep_jumps = [[1]]
jump_to_dep_specs = [[1]]
dg = [[1]]
majump = MassActionJump(rates, reactstoch, netstoch)
dprob = DiscreteProblem(u0,(0.,10.),rates)

algs = DiffEqJump.JUMP_AGGREGATORS

for ssa in algs
    jprob = JumpProblem(dprob, ssa, majump; vartojumps_map=spec_to_dep_jumps, jumptovars_map=jump_to_dep_specs)
    sol = solve(jprob, SSAStepper())
    @test sol[1,end] == 0
    @test sol.t[end] < Inf
end


text:
using jcontrol
import jcontrol.comparisons: track_segments
using jcontrol.visuals
import jcontrol: trimtrial
using Plots

"""
Load tracking data and cut between start and end of selected
segment, then visualize tracking and other variables.
"""

segment = track_segments[3]

# load & trim trials
trials = @time load_cached_trials(; keep_n=20)
trials = map(t->trimtrial(t, segment), trials)

varvalues = Dict()
for var in (:x, :y, :θ, :u :ω)
    varvalues[var] = Dict(
        "start" => get_varvalue_at_frameidx(trials, var, 1),
        "stop" => get_varvalue_at_frameidx(trials, var)
    )
end

# TODO draw stuff
# ------------------------------- visualization ------------------------------ #
# draw tracking data
p1 = draw(:arena)
draw!(track; alpha=.1)
draw!(segment)
draw!.(trials)

function plotdist(key)
    hist = histogram(varvalues[key]["start"], label=string(key)*" init")
    hist = histogram!(varvalues[key]["stop"], label=string(key)*" final")
    return hist
end

plots = map(plotdist, collect(keys(varvalues)))

l = @layout [
    a{0.5w} grid(length(keys(varvalues)), 1)
]

closeall()
display(
    plot(
        p1,
        plots...;
        layout=l,
    )
)

max_stars_repo_name: FedeClaudi/LocomotionControl

max_issues_repo_name: FedeClaudi/LocomotionControl

max_forks_repo_name: FedeClaudi/LocomotionControl


text:
using Test
using ADBCVUMPS
using ADBCVUMPS:num_grad
using BCVUMPS
using BCVUMPS:bcvumps_env,magnetisation,Z
using Random
using Zygote

@testset "$(Ni)x$(Nj) ising" for Ni = [1,2], Nj = [1,2]
    Random.seed!(100)
    model = Ising(Ni,Nj)
    β,D = 0.5, 10
    function foo1(β)
        -log(Z(bcvumps_env(model,β,D)))
    end
    function foo2(β)
        magnetisation(bcvumps_env(model,β,D; maxiter = 20, tol = 1e-12),model,β)
    end

    @test isapprox(Zygote.gradient(foo1,β)[1], num_grad(foo1,β), atol = 1e-6)
    @test isapprox(Zygote.gradient(foo2,β)[1], num_grad(foo2,β), atol = 1e-6)
end

@testset "J1-J2-2x2-ising" begin
    Random.seed!(100)
    model = Ising22(2.0)
    β,D = 0.5,10
    function foo1(β)
        -log(Z(bcvumps_env(model,β,D)))
    end
    @test isapprox(Zygote.gradient(foo1,β)[1], num_grad(foo1,β), atol = 1e-6)
end


text:
export NNet, evaluate_network, evaluate_network_multiple, num_inputs, num_outputs

"""
Custom type that represents a fully connected ReLU network from a .nnet file

Args:
    file (string): A .nnet file to load

Attributes:
    numLayers (int): Number of weight matrices or bias vectors in neural network
    layerSizes (list of ints): Size of input layer, hidden layers, and output layer
    inputSize (int): Size of input
    outputSize (int): Size of output
    mins (list of floats): Minimum values of inputs
    maxes (list of floats): Maximum values of inputs
    means (list of floats): Means of inputs and mean of outputs
    ranges (list of floats): Ranges of inputs and range of outputs
    weights (list of arrays): Weight matrices in network
    biases (list of arrays): Bias vectors in network
"""
type NNet
    weights::Array{Any,1}
    biases::Array{Any,1}
    numLayers::Int32
    layerSizes::Array{Int32,1}
    inputSize::Int32
    outputSize::Int32
    mins::Array{Float64,1}
    maxes::Array{Float64,1}
    means::Array{Float64,1}
    ranges::Array{Float64,1}

    function NNet(file::AbstractString)
        this  = new()

        # Open file for reading
        f = open(file)

        # Skip any header lines
        line = readline(f)
        while line[1:2]=="//"
            line = readline(f)
        end

        # Read information about the neural network
        record = split(line,[',','\n'])
        this.numLayers = parse(Int32,record[1])
        this.inputSize = parse(Int32,record[2])
        this.outputSize = parse(Int32,record[3])

        line = readline(f)
        record = split(line,[',','\n'])
        this.layerSizes = zeros(this.numLayers+1)
        for i=1:(this.numLayers+1)
            this.layerSizes[i]=parse(Int32,record[i])
        end

        # Skip unused line
        line = readline(f)

        line = readline(f)
        record = split(line,[',','\n'])
        this.mins = zeros(this.inputSize)
        for i=1:(this.inputSize)
            this.mins[i]=parse(Float64,record[i])
        end

        line = readline(f)
        record = split(line,[',','\n'])
        this.maxes = zeros(this.inputSize)
        for i=1:(this.inputSize)
            this.maxes[i]=parse(Float64,record[i])
        end

        line = readline(f)
        record = split(line,[',','\n'])
        this.means = zeros(this.inputSize+1)
        for i=1:(this.inputSize+1)
            this.means[i]=parse(Float64,record[i])
        end

        line = readline(f)
        record = split(line,[',','\n'])
        this.ranges = zeros(this.inputSize+1)
        for i=1:(this.inputSize+1)
            this.ranges[i]=parse(Float64,record[i])
        end

        # Initialize weight and bias arrays
        this.weights = Any[zeros(this.layerSizes[2],this.layerSizes[1])]
        this.biases  = Any[zeros(this.layerSizes[2])]
        for i=2:this.numLayers
            this.weights = [this.weights;Any[zeros(this.layerSizes[i+1],this.layerSizes[i])]]
            this.biases  = [this.biases;Any[zeros(this.layerSizes[i+1])]]
        end

        # Fill weight and bias arrays with values from nnet file
        layer=1
        i=1
        j=1
        line = readline(f)
        record = split(line,[',','\n'])
        while !eof(f)
            while i<=this.layerSizes[layer+1]
                while record[j]!=""
                    this.weights[layer][i,j] = parse(Float64,record[j])
                    j=j+1
                end
                j=1
                i=i+1
                line = readline(f)
                record = split(line,[',','\n'])
            end
            i=1
            while i<=this.layerSizes[layer+1]
                this.biases[layer][i] = parse(Float64,record[1])
                i=i+1
                line = readline(f)
                record = split(line,[',','\n'])
            end
            layer=layer+1
            i=1
            j=1
        end
        close(f)

        return this
    end
end

"""
Evaluate network using given inputs

Args:
    nnet (NNet): Neural network to evaluate
    inputs (array): Network inputs to be evaluated

Returns:
    (array): Network output
"""
function evaluate_network(nnet::NNet,input::Array{Float64,1})
    numLayers = nnet.numLayers
    inputSize = nnet.inputSize
    outputSize = nnet.outputSize
    biases = nnet.biases
    weights = nnet.weights

    # Prepare the inputs to the neural network
    inputs = zeros(inputSize)
    for i = 1:inputSize
        if input[i]<nnet.mins[i]
            inputs[i] = (nnet.mins[i]-nnet.means[i])/nnet.ranges[i]
        elseif input[i] > nnet.maxes[i]
            inputs[i] = (nnet.maxes[i]-nnet.means[i])/nnet.ranges[i]
        else
            inputs[i] = (input[i]-nnet.means[i])/nnet.ranges[i]
        end
    end

    # Evaluate the neural network
    for layer = 1:numLayers-1
        temp = max.(*(weights[layer],inputs[1:nnet.layerSizes[layer]])+biases[layer],0)
        inputs = temp
    end
    outputs = *(weights[end],inputs[1:nnet.layerSizes[end-1]])+biases[end]

    # Undo output normalization
    for i=1:outputSize
        outputs[i] = outputs[i]*nnet.ranges[end]+nnet.means[end]
    end
    return outputs
end

"""
Evaluate network using multiple sets of inputs

Args:
    nnet (NNet): Neural network to evaluate
    inputs (array): Network inputs to be evaluated

Returns:
    (array): Network outputs for each set of inputs
"""
function evaluate_network_multiple(nnet::NNet,input::Array{Float64,2})
    numLayers = nnet.numLayers
    inputSize = nnet.inputSize
    outputSize = nnet.outputSize
    biases = nnet.biases
    weights = nnet.weights

    # Prepare the inputs to the neural network
    _,numInputs = size(input)
    inputs = zeros(inputSize,numInputs)
    for i = 1:inputSize
        for j = 1:numInputs
            if input[i,j]<nnet.mins[i]
                inputs[i,j] = (nnet.mins[i]-nnet.means[i])/nnet.ranges[i]
            elseif input[i,j] > nnet.maxes[i]
                inputs[i,j] = (nnet.maxes[i]-nnet.means[i])/nnet.ranges[i]
            else
                inputs[i,j] = (input[i,j]-nnet.means[i])/nnet.ranges[i]
            end
        end
    end

    # Evaluate the neural network
    for layer = 1:numLayers-1
        inputs = max.(*(weights[layer],inputs[1:nnet.layerSizes[layer],:])+*(biases[layer],ones(1,numInputs)),0)
    end
    outputs = *(weights[end],inputs[1:nnet.layerSizes[end-1],:])+*(biases[end],ones(1,numInputs))

    # Undo output normalization
    for i=1:outputSize
        for j=1:numInputs
            outputs[i,j] = outputs[i,j]*nnet.ranges[end]+nnet.means[end]
        end
    end
    return outputs
end

""" Get number of inputs to network"""
function num_inputs(nnet::NNet)
    return nnet.inputSize
end

""" Get number of outputs from network"""
function num_outputs(nnet::NNet)
    return nnet.outputSize
end


text:
"""
    AbstractIndexOrder

Abstract type determining the ordering of dimensions of an array with possibly
permuted indices.

Subtypes are [`MemoryOrder`](@ref) and [`LogicalOrder`](@ref).
"""
abstract type AbstractIndexOrder end

"""
    MemoryOrder <: AbstractIndexOrder

Singleton type specifying that array dimensions should be given in memory (or
*permuted*) order.
"""
struct MemoryOrder <: AbstractIndexOrder end

"""
    LogicalOrder <: AbstractIndexOrder

Singleton type specifying that array dimensions should be given in logical (or
*non-permuted*) order.
"""
struct LogicalOrder <: AbstractIndexOrder end

const DefaultOrder = LogicalOrder


text:
struct GridData
    voxelNums::Vector{Int}
    startVox::Vector{Float64}
    voxelWidths::Vector{Float64}
    endVox::Vector{Float64}
end

struct ROIData
    Mins::Vector{Float64}
    Maxs::Vector{Float64}
end

function createGrid( minPoint, maxPoint, voxelNums )

    @assert sum( minPoint .<= maxPoint ) == 3

    ### MODIFY USING MULTIPLE DISPATCH ###
    if length( voxelNums ) == 1
        voxelNums = voxelNums * ones(3)
    end
    # ----------------------------------

    voxelWidths   = ( maxPoint - minPoint ) ./ voxelNums
    startVox      =   minPoint + voxelWidths * 0.5
    endVox        =   maxPoint - voxelWidths * 0.5

    return GridData( voxelNums, startVox, voxelWidths, endVox )
end


text:
@testset "Simple generation and general methods" begin
  s = CD.Sample(CD.Size(23,24), 23, [0, 1])
  @test length(CD.patches(s)) == 0
  @test size(s) == (23,24)
  @test size(s,1) == 23
  @test size(s,2) == 24
end

@testset "Random Patch fill" begin
  s = CD.Sample(CD.Size(64), 2, [0, 1])
  s = rand(s; patchsize=CD.Size(8), nbpatches=5, nbsame=3, diff=2)
  @test length(CD.patches(s)) == 5
  for i in CD.patches(s)
    @test 0 ∈ CD.content(i)
    @test 1 ∈ CD.content(i)
    @test 2 ∉ CD.content(i)
  end

  rd = CD.render(s)
  @test 0 ∈ rd
  @test 1 ∈ rd
  @test 2 ∈ rd
  @test size(rd) == (64, 64)
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# Methods operating on different special matrix types

# Interconversion between special matrix types

# conversions from Diagonal to other special matrix types
Bidiagonal(A::Diagonal) = Bidiagonal(A.diag, fill!(similar(A.diag, length(A.diag)-1), 0), :U)
SymTridiagonal(A::Diagonal) = SymTridiagonal(A.diag, fill!(similar(A.diag, length(A.diag)-1), 0))
Tridiagonal(A::Diagonal) = Tridiagonal(fill!(similar(A.diag, length(A.diag)-1), 0), A.diag,
                                       fill!(similar(A.diag, length(A.diag)-1), 0))

# conversions from Bidiagonal to other special matrix types
Diagonal(A::Bidiagonal) =
    iszero(A.ev) ? Diagonal(A.dv) :
        throw(ArgumentError("matrix cannot be represented as Diagonal"))
SymTridiagonal(A::Bidiagonal) =
    iszero(A.ev) ? SymTridiagonal(A.dv, A.ev) :
        throw(ArgumentError("matrix cannot be represented as SymTridiagonal"))
Tridiagonal(A::Bidiagonal) =
    Tridiagonal(A.uplo == 'U' ? fill!(similar(A.ev), 0) : A.ev, A.dv,
                A.uplo == 'U' ? A.ev : fill!(similar(A.ev), 0))

# conversions from SymTridiagonal to other special matrix types
Diagonal(A::SymTridiagonal) =
    iszero(A.ev) ? Diagonal(A.dv) :
        throw(ArgumentError("matrix cannot be represented as Diagonal"))
Bidiagonal(A::SymTridiagonal) =
    iszero(A.ev) ? Bidiagonal(A.dv, A.ev, :U) :
        throw(ArgumentError("matrix cannot be represented as Bidiagonal"))
Tridiagonal(A::SymTridiagonal) =
    Tridiagonal(copy(A.ev), A.dv, A.ev)

# conversions from Tridiagonal to other special matrix types
Diagonal(A::Tridiagonal) =
    iszero(A.dl) && iszero(A.du) ? Diagonal(A.d) :
        throw(ArgumentError("matrix cannot be represented as Diagonal"))
Bidiagonal(A::Tridiagonal) =
    iszero(A.dl) ? Bidiagonal(A.d, A.du, :U) :
    iszero(A.du) ? Bidiagonal(A.d, A.dl, :L) :
        throw(ArgumentError("matrix cannot be represented as Bidiagonal"))
SymTridiagonal(A::Tridiagonal) =
    A.dl == A.du ? SymTridiagonal(A.d, A.dl) :
        throw(ArgumentError("matrix cannot be represented as SymTridiagonal"))

# conversions from AbstractTriangular to special matrix types
Diagonal(A::AbstractTriangular) =
    isdiag(A) ? Diagonal(diag(A)) :
        throw(ArgumentError("matrix cannot be represented as Diagonal"))
Bidiagonal(A::AbstractTriangular) =
    isbanded(A, 0, 1) ? Bidiagonal(diag(A, 0), diag(A,  1), :U) : # is upper bidiagonal
    isbanded(A, -1, 0) ? Bidiagonal(diag(A, 0), diag(A, -1), :L) : # is lower bidiagonal
        throw(ArgumentError("matrix cannot be represented as Bidiagonal"))
SymTridiagonal(A::AbstractTriangular) = SymTridiagonal(Tridiagonal(A))
Tridiagonal(A::AbstractTriangular) =
    isbanded(A, -1, 1) ? Tridiagonal(diag(A, -1), diag(A, 0), diag(A, 1)) : # is tridiagonal
        throw(ArgumentError("matrix cannot be represented as Tridiagonal"))

const ConvertibleSpecialMatrix = Union{Diagonal,Bidiagonal,SymTridiagonal,Tridiagonal,AbstractTriangular}

convert(T::Type{<:Diagonal},       m::ConvertibleSpecialMatrix) = m isa T ? m : T(m)
convert(T::Type{<:SymTridiagonal}, m::ConvertibleSpecialMatrix) = m isa T ? m : T(m)
convert(T::Type{<:Tridiagonal},    m::ConvertibleSpecialMatrix) = m isa T ? m : T(m)

convert(T::Type{<:LowerTriangular}, m::Union{LowerTriangular,UnitLowerTriangular}) = m isa T ? m : T(m)
convert(T::Type{<:UpperTriangular}, m::Union{UpperTriangular,UnitUpperTriangular}) = m isa T ? m : T(m)

# Constructs two method definitions taking into account (assumed) commutativity
# e.g. @commutative f(x::S, y::T) where {S,T} = x+y is the same is defining
#     f(x::S, y::T) where {S,T} = x+y
#     f(y::T, x::S) where {S,T} = f(x, y)
macro commutative(myexpr)
    @assert myexpr.head===:(=) || myexpr.head===:function # Make sure it is a function definition
    y = copy(myexpr.args[1].args[2:end])
    reverse!(y)
    reversed_call = Expr(:(=), Expr(:call,myexpr.args[1].args[1],y...), myexpr.args[1])
    esc(Expr(:block, myexpr, reversed_call))
end

for op in (:+, :-)
    SpecialMatrices = [:Diagonal, :Bidiagonal, :Tridiagonal, :Matrix]
    for (idx, matrixtype1) in enumerate(SpecialMatrices) # matrixtype1 is the sparser matrix type
        for matrixtype2 in SpecialMatrices[idx+1:end] # matrixtype2 is the denser matrix type
            @eval begin # TODO quite a few of these conversions are NOT defined
                ($op)(A::($matrixtype1), B::($matrixtype2)) = ($op)(convert(($matrixtype2), A), B)
                ($op)(A::($matrixtype2), B::($matrixtype1)) = ($op)(A, convert(($matrixtype2), B))
            end
        end
    end

    for  matrixtype1 in (:SymTridiagonal,)                      # matrixtype1 is the sparser matrix type
        for matrixtype2 in (:Tridiagonal, :Matrix) # matrixtype2 is the denser matrix type
            @eval begin
                ($op)(A::($matrixtype1), B::($matrixtype2)) = ($op)(convert(($matrixtype2), A), B)
                ($op)(A::($matrixtype2), B::($matrixtype1)) = ($op)(A, convert(($matrixtype2), B))
            end
        end
    end

    for matrixtype1 in (:Diagonal, :Bidiagonal) # matrixtype1 is the sparser matrix type
        for matrixtype2 in (:SymTridiagonal,)   # matrixtype2 is the denser matrix type
            @eval begin
                ($op)(A::($matrixtype1), B::($matrixtype2)) = ($op)(convert(($matrixtype2), A), B)
                ($op)(A::($matrixtype2), B::($matrixtype1)) = ($op)(A, convert(($matrixtype2), B))
            end
        end
    end

    for matrixtype1 in (:Diagonal,)
        for (matrixtype2,matrixtype3) in ((:UpperTriangular,:UpperTriangular),
                                          (:UnitUpperTriangular,:UpperTriangular),
                                          (:LowerTriangular,:LowerTriangular),
                                          (:UnitLowerTriangular,:LowerTriangular))
            @eval begin
                ($op)(A::($matrixtype1), B::($matrixtype2)) = ($op)(($matrixtype3)(A), B)
                ($op)(A::($matrixtype2), B::($matrixtype1)) = ($op)(A, ($matrixtype3)(B))
            end
        end
    end
    for matrixtype in (:SymTridiagonal,:Tridiagonal,:Bidiagonal,:Matrix)
        @eval begin
            ($op)(A::AbstractTriangular, B::($matrixtype)) = ($op)(copyto!(similar(parent(A)), A), B)
            ($op)(A::($matrixtype), B::AbstractTriangular) = ($op)(A, copyto!(similar(parent(B)), B))
        end
    end
end

rmul!(A::AbstractTriangular, adjB::Adjoint{<:Any,<:Union{QRCompactWYQ,QRPackedQ}}) =
    (B = adjB.parent; rmul!(full!(A), adjoint(B)))
*(A::AbstractTriangular, adjB::Adjoint{<:Any,<:Union{QRCompactWYQ,QRPackedQ}}) =
    (B = adjB.parent; *(copyto!(similar(parent(A)), A), adjoint(B)))

# fill[stored]! methods
fillstored!(A::Diagonal, x) = (fill!(A.diag, x); A)
fillstored!(A::Bidiagonal, x) = (fill!(A.dv, x); fill!(A.ev, x); A)
fillstored!(A::Tridiagonal, x) = (fill!(A.dl, x); fill!(A.d, x); fill!(A.du, x); A)
fillstored!(A::SymTridiagonal, x) = (fill!(A.dv, x); fill!(A.ev, x); A)

_small_enough(A::Union{Diagonal, Bidiagonal}) = size(A, 1) <= 1
_small_enough(A::Tridiagonal) = size(A, 1) <= 2
_small_enough(A::SymTridiagonal) = size(A, 1) <= 2

function fill!(A::Union{Diagonal,Bidiagonal,Tridiagonal,SymTridiagonal}, x)
    xT = convert(eltype(A), x)
    (iszero(xT) || _small_enough(A)) && return fillstored!(A, xT)
    throw(ArgumentError("array of type $(typeof(A)) and size $(size(A)) can
    not be filled with $x, since some of its entries are constrained."))
end

# equals and approx equals methods for structured matrices
# SymTridiagonal == Tridiagonal is already defined in tridiag.jl

# SymTridiagonal and Bidiagonal have the same field names
==(A::Diagonal, B::Union{SymTridiagonal, Bidiagonal}) = iszero(B.ev) && A.diag == B.dv
==(B::Bidiagonal, A::Diagonal) = A == B

==(A::Diagonal, B::Tridiagonal) = iszero(B.dl) && iszero(B.du) && A.diag == B.d
==(B::Tridiagonal, A::Diagonal) = A == B

function ==(A::Bidiagonal, B::Tridiagonal)
    if A.uplo == 'U'
        return iszero(B.dl) && A.dv == B.d && A.ev == B.du
    else
        return iszero(B.du) && A.dv == B.d && A.ev == B.dl
    end
end
==(B::Tridiagonal, A::Bidiagonal) = A == B

==(A::Bidiagonal, B::SymTridiagonal) = iszero(B.ev) && iszero(A.ev) && A.dv == B.dv
==(B::SymTridiagonal, A::Bidiagonal) = A == B

max_stars_repo_name: mhaldrich/julia

max_issues_repo_name: mhaldrich/julia

max_forks_repo_name: mhaldrich/julia


text:
# This file was generated, do not modify it. # hide
best_model = fitted_params(mtm).best_model
@show best_model.learning_rate
@show best_model.min_data_in_leaf
@show best_model.num_iterations


text:
using Clapeyron, Test, NLSolvers, ForwardDiff

const SOL = Clapeyron.Solvers

quadratic(x) = x*x - 4
rosenbrock(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
golden_number_fixpoint(x) = one(x) + one(x)/x
fgh_lgmx(x) = (log(x)+x,1/x +1,-1/(x^2))
function quadratic_fixpoint(y,x)
    y[1] = sqrt(-x[2]^2 + 5)
    y[2] = 2*sqrt(-x[1]^2 + 9)/3
    y
end
@testset "Solvers Module" begin
    @testset "newton,halley" begin
        x0 = 2+rand()
        fdf = SOL.f∂f(quadratic,x0)
        @test fdf[1] ≈ quadratic(x0)
        @test fdf[2] ≈ 2*x0
        @test SOL.ad_newton(quadratic,x0) ≈ 2.0
        @test SOL.newton(x->(quadratic(x),2*x),x0) ≈ 2.0
        @test SOL.halley(fgh_lgmx,0.5)≈ 0.567143290409784
    end

    @testset "box optimize" begin
        #example of FMinBox in Optim.jl
        lower = [1.25, -2.1]
        upper = [Inf, Inf]
        solution = [1.25, 1.5625]
        initial_x = [2.0, 2.0]
        res = SOL.box_optimize(rosenbrock,initial_x,lower,upper)
        @test all(res.info.minimizer .≈ solution)
        @test res.info.minimum ≈ 0.0625
    end

    @testset "fixpoint" begin
        #example of FMinBox in Optim.jl
        x0 = 2+ rand()
        ϕ = (1+sqrt(5))/2
        x1 = [0.1,0.1]
        @test @inferred(SOL.fixpoint(golden_number_fixpoint,x0)) ≈ ϕ
        @test @inferred(SOL.fixpoint(golden_number_fixpoint,x0,SOL.AitkenFixPoint())) ≈ ϕ
        @test @inferred(SOL.fixpoint(quadratic_fixpoint,x1)) ≈ [0.6*sqrt(5),0.8*sqrt(5)]
    end

    function f_diffmcp!(fvec, x)
        fvec[1] = (1-x[1])^2-1.01
    end

    @testset "roots3" begin
        #@test [@inferred(SOL.roots3(SA[-6im, -(3 + 4im), 2im-2, 1.0]))...] ≈ [3, -2im, -1]
        @test @inferred(SOL.roots3(1.0, -3.0, 3.0, -1.0)) ≈ [1, 1, 1]
    end
    # A difficult MCP.
    #
    # Presented in Miranda and Fackler (2002): "Applied Computational Economics and
    # Finance", p. 51
    #does not converge in NLSolve.jl converges here with NLSolvers.jl
    @testset "nlsolve" begin
        res = SOL.nlsolve(f_diffmcp!,[0.0],TrustRegion(Newton()),ForwardDiff.Chunk{1}())
        res2 = SOL.nlsolve(f_diffmcp!,[0.0],LineSearch(Newton()),ForwardDiff.Chunk{1}())
        @test SOL.x_sol(res) isa Vector
        @test SOL.x_sol(res2) isa Vector
        solution = SOL.x_sol(res)
        v = [1.0]
        f_diffmcp!(v,solution)
        @test v[1] <= sqrt(eps(Float64))
    end

    @testset "Rachford Rice" begin
    #error, all Ks > 0, from CalebBell/Chemicals
    zs = [0.0, 0.0885053990596404, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.721469918219507, 0.01742948033685831,
            0.1725952023839942]
    Ks = [192.3625321718047, 105.20070698573475, 76.30532397111489, 37.090890262982, 21.38862102676539, 18.093547012968767,
            10.319129068837443, 9.001962137200403, 4.565198737490148, 340.69153314749224, 269.09234343328467,
            21.858385052861507]
    @test Clapeyron.rr_vle_vapor_fraction(Ks,zs) == Inf

    #here are some tests, from the paper.
    #the paper has errors, z2 has 10 components, 9 and 10 repeated.
    #and after that, they dont include the bmax,bmin in their selected pair.
    z1 = [0.30,0.15, 0.05, 0.02, 0.01, 0.02, 0.02, 0.03, 0.07, 0.33]
    k1 = [3.0E+00, 2.0E+00, 1.1E+00, 8.0E-01, 4.0E-01, 1.0E-01, 5.0E-02, 2.0E-02, 1.0E-02, 1.0E-04]
    idxs1,ks1,zs1 = Clapeyron.rr_find_strongest(k1,z1)
    @test all(in.(idxs1,Ref((1,2,9,10))))
    @test Clapeyron.rr_vle_vapor_fraction(k1,z1) ≈ 0.1769 atol = 1e-4
    z2 = [0.00034825,0.01376300,0.13084000,0.10925000,0.00001000,0.51009000,0.23564000,0.00006000]
    k2 = [5.2663E+02, 5.0400E+01, 1.6463E+00, 8.7450E-01, 1.5589E-01, 3.6588E-02, 2.6625E-02, 4.8918E-06]
    idxs2,ks2,zs2 = Clapeyron.rr_find_strongest(k2,z2)
    #the paper is wrong here, it stablishes (1,2,6,7) ,does not include bmax
    #we still converge to the correct root
    @test all(in.(idxs2,Ref((1,2,6,8))))
    @test Clapeyron.rr_vle_vapor_fraction(k2,z2) ≈ 0.00327 atol = 1e-4

    z3 = [0.0187002, 0.0243002, 0.5419054, 0.0999010, 0.0969010, 0.0400004, 0.0212002, 0.0148001, 0.0741507, 0.0350404, 0.0173602, 0.0157402]
    k3 = [1.32420, 1.12778, 1.22222, 1.11760, 9.88047E-01, 8.94344E-01, 7.87440E-01, 7.43687E-01, 8.11797E-01, 6.93279E-01, 5.09443E-01, 2.28721E-01]
    idxs3,ks3,zs3 = Clapeyron.rr_find_strongest(k3,z3)
    @test all(in.(idxs3,Ref((1,3,9,12))))
    @test Clapeyron.rr_vle_vapor_fraction(k3,z3) ≈ 0.98725 atol = 1e-4

    #testing exact solver,3 comps:
    zs_3 = [0.5, 0.3, 0.2]
    Ks_3 = [1.685, 0.742, 0.532]
    xs_3 = [0.33940869696634357, 0.3650560590371706, 0.2955352439964858]
    ys_3 = [0.5719036543882889, 0.27087159580558057, 0.15722474980613044]
    β_3 = 0.6907302627738544
    β_3_sol = Clapeyron.rr_vle_vapor_fraction(Ks_3,zs_3)
    @test  β_3_sol ≈ β_3
    @test Clapeyron.rr_flash_eval(Ks_3,zs_3,β_3_sol) <= 4*eps(β_3)
    @test Clapeyron.rr_flash_vapor(Ks_3,zs_3,β_3_sol) ≈ ys_3
    @test Clapeyron.rr_flash_liquid(Ks_3,zs_3,β_3_sol) ≈ xs_3

    #Extreme Kvalues, it should give β = 0.99999
    Ks_extreme = [15.464909530837806, 7006.64008090944, 1.8085837711444488, 0.007750676421035811, 30.98450366497431]
    zs_extreme = [0.26562380186293233, 0.04910234829974003, 0.284394553603828, 0.006300023876072552, 0.3945792723574272]
    @test Clapeyron.rr_vle_vapor_fraction(Ks_extreme,zs_extreme) ≈ 0.999999999

    # Case where the evaluated point is right on the boundary
    Ks_eps_0 = [1.2566703532018493e-21, 3.35062752053393, 1.0300675710905643e-23, 1.706258568414198e-39, 1.6382855298440747e-20]
    zs_eps_0 = [0.13754371891028325, 0.2984515568715462, 0.2546683930289046, 0.08177453852283137, 0.22756179266643456]
    @test abs(Clapeyron.rr_vle_vapor_fraction(Ks_eps_0,zs_eps_0)) < eps(Float64)
end
    @testset "det_22" begin
        a1 = 1
        a2 = 2
        a3 = 3
        a4 = 4
        @test SOL.det_22(1,2,3,4) == a1*a2 - a3*a4
    end
end


text:
using BenchmarkTools
using LightGraphs

DIGRAPHS = Dict{String,DiGraph}(
    "complete100"   => complete_digraph(100),
    "path500"       => path_digraph(500)
)

GRAPHS = Dict{String,Graph}(
    "complete100"   => complete_graph(100),
    "tutte"         => smallgraph(:tutte),
    "path500"       => path_graph(500)
)

suite = BenchmarkGroup()
include("core.jl")

tune!(suite);
results = run(suite, verbose = true, seconds = 10)