text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -37.14009556780016 <= q <= 48.365016418887166)

                     add_NL_constraint(m, :(log(1 + exp(0.24902435571631587 + -0.02312729681609138*log(1 + exp(-0.46870412880299783 + 0.1700038462852933*log(1 + exp(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3]))) + -0.44738483682757924*log(1 + exp(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3]))) + -0.5112184633122538*log(1 + exp(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3]))) + 0.17714072857386487*log(1 + exp(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3]))) + 0.9529511507705823*log(1 + exp(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))))) + -0.5418072450548097*log(1 + exp(0.5019192867350393 + 0.2997113620022893*log(1 + exp(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3]))) + -0.5721484689942011*log(1 + exp(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3]))) + 0.4073624372378828*log(1 + exp(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3]))) + -0.07630883358086304*log(1 + exp(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3]))) + 0.3103910161241421*log(1 + exp(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))))) + 0.9789829058919253*log(1 + exp(-0.652369829242545 + -0.998242239199818*log(1 + exp(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3]))) + -0.041041478696029365*log(1 + exp(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3]))) + -0.4348093907280841*log(1 + exp(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3]))) + 0.7773011179238121*log(1 + exp(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3]))) + 0.7537599741964742*log(1 + exp(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))))) + -0.13446361488468206*log(1 + exp(-0.40141742757451926 + 0.7289092837988314*log(1 + exp(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3]))) + 0.16077592157977172*log(1 + exp(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3]))) + 0.6378484913808045*log(1 + exp(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3]))) + 0.9326240829024273*log(1 + exp(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3]))) + 0.9638370498324353*log(1 + exp(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))))) + -0.3069132295129271*log(1 + exp(0.17680227033398044 + 0.8293506390949057*log(1 + exp(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3]))) + -0.09356004658531791*log(1 + exp(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3]))) + -0.16358082594109513*log(1 + exp(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3]))) + 0.9029174038261001*log(1 + exp(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3]))) + 0.7318887196209252*log(1 + exp(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))))))) + log(1 + exp(-0.017552134571214406 + -0.41004797399293524*log(1 + exp(-0.46870412880299783 + 0.1700038462852933*log(1 + exp(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3]))) + -0.44738483682757924*log(1 + exp(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3]))) + -0.5112184633122538*log(1 + exp(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3]))) + 0.17714072857386487*log(1 + exp(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3]))) + 0.9529511507705823*log(1 + exp(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))))) + -0.47019448283071075*log(1 + exp(0.5019192867350393 + 0.2997113620022893*log(1 + exp(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3]))) + -0.5721484689942011*log(1 + exp(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3]))) + 0.4073624372378828*log(1 + exp(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3]))) + -0.07630883358086304*log(1 + exp(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3]))) + 0.3103910161241421*log(1 + exp(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))))) + 0.35795867979675133*log(1 + exp(-0.652369829242545 + -0.998242239199818*log(1 + exp(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3]))) + -0.041041478696029365*log(1 + exp(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3]))) + -0.4348093907280841*log(1 + exp(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3]))) + 0.7773011179238121*log(1 + exp(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3]))) + 0.7537599741964742*log(1 + exp(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))))) + 0.44715886126593407*log(1 + exp(-0.40141742757451926 + 0.7289092837988314*log(1 + exp(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3]))) + 0.16077592157977172*log(1 + exp(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3]))) + 0.6378484913808045*log(1 + exp(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3]))) + 0.9326240829024273*log(1 + exp(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3]))) + 0.9638370498324353*log(1 + exp(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))))) + -0.22638221307132733*log(1 + exp(0.17680227033398044 + 0.8293506390949057*log(1 + exp(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3]))) + -0.09356004658531791*log(1 + exp(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3]))) + -0.16358082594109513*log(1 + exp(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3]))) + 0.9029174038261001*log(1 + exp(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3]))) + 0.7318887196209252*log(1 + exp(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))))))) + log(1 + exp(-0.4029592795849051 + -0.6197606146665708*log(1 + exp(-0.46870412880299783 + 0.1700038462852933*log(1 + exp(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3]))) + -0.44738483682757924*log(1 + exp(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3]))) + -0.5112184633122538*log(1 + exp(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3]))) + 0.17714072857386487*log(1 + exp(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3]))) + 0.9529511507705823*log(1 + exp(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))))) + 0.8464388391219866*log(1 + exp(0.5019192867350393 + 0.2997113620022893*log(1 + exp(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3]))) + -0.5721484689942011*log(1 + exp(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3]))) + 0.4073624372378828*log(1 + exp(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3]))) + -0.07630883358086304*log(1 + exp(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3]))) + 0.3103910161241421*log(1 + exp(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))))) + 0.6262663939568953*log(1 + exp(-0.652369829242545 + -0.998242239199818*log(1 + exp(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3]))) + -0.041041478696029365*log(1 + exp(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3]))) + -0.4348093907280841*log(1 + exp(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3]))) + 0.7773011179238121*log(1 + exp(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3]))) + 0.7537599741964742*log(1 + exp(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))))) + -0.06090941414894013*log(1 + exp(-0.40141742757451926 + 0.7289092837988314*log(1 + exp(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3]))) + 0.16077592157977172*log(1 + exp(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3]))) + 0.6378484913808045*log(1 + exp(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3]))) + 0.9326240829024273*log(1 + exp(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3]))) + 0.9638370498324353*log(1 + exp(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))))) + -0.6699649674430836*log(1 + exp(0.17680227033398044 + 0.8293506390949057*log(1 + exp(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3]))) + -0.09356004658531791*log(1 + exp(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3]))) + -0.16358082594109513*log(1 + exp(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3]))) + 0.9029174038261001*log(1 + exp(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3]))) + 0.7318887196209252*log(1 + exp(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))))))) + log(1 + exp(0.6579858782021573 + -0.21911632564812766*log(1 + exp(-0.46870412880299783 + 0.1700038462852933*log(1 + exp(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3]))) + -0.44738483682757924*log(1 + exp(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3]))) + -0.5112184633122538*log(1 + exp(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3]))) + 0.17714072857386487*log(1 + exp(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3]))) + 0.9529511507705823*log(1 + exp(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))))) + 0.7292607357990395*log(1 + exp(0.5019192867350393 + 0.2997113620022893*log(1 + exp(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3]))) + -0.5721484689942011*log(1 + exp(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3]))) + 0.4073624372378828*log(1 + exp(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3]))) + -0.07630883358086304*log(1 + exp(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3]))) + 0.3103910161241421*log(1 + exp(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))))) + -0.23882038501163416*log(1 + exp(-0.652369829242545 + -0.998242239199818*log(1 + exp(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3]))) + -0.041041478696029365*log(1 + exp(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3]))) + -0.4348093907280841*log(1 + exp(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3]))) + 0.7773011179238121*log(1 + exp(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3]))) + 0.7537599741964742*log(1 + exp(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))))) + 0.6855544316424598*log(1 + exp(-0.40141742757451926 + 0.7289092837988314*log(1 + exp(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3]))) + 0.16077592157977172*log(1 + exp(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3]))) + 0.6378484913808045*log(1 + exp(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3]))) + 0.9326240829024273*log(1 + exp(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3]))) + 0.9638370498324353*log(1 + exp(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))))) + -0.2868272308560216*log(1 + exp(0.17680227033398044 + 0.8293506390949057*log(1 + exp(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3]))) + -0.09356004658531791*log(1 + exp(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3]))) + -0.16358082594109513*log(1 + exp(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3]))) + 0.9029174038261001*log(1 + exp(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3]))) + 0.7318887196209252*log(1 + exp(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))))))) + log(1 + exp(-0.10147744083820065 + -0.7175110849623545*log(1 + exp(-0.46870412880299783 + 0.1700038462852933*log(1 + exp(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3]))) + -0.44738483682757924*log(1 + exp(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3]))) + -0.5112184633122538*log(1 + exp(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3]))) + 0.17714072857386487*log(1 + exp(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3]))) + 0.9529511507705823*log(1 + exp(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))))) + 0.9583278917345623*log(1 + exp(0.5019192867350393 + 0.2997113620022893*log(1 + exp(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3]))) + -0.5721484689942011*log(1 + exp(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3]))) + 0.4073624372378828*log(1 + exp(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3]))) + -0.07630883358086304*log(1 + exp(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3]))) + 0.3103910161241421*log(1 + exp(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))))) + 0.7313133539296679*log(1 + exp(-0.652369829242545 + -0.998242239199818*log(1 + exp(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3]))) + -0.041041478696029365*log(1 + exp(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3]))) + -0.4348093907280841*log(1 + exp(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3]))) + 0.7773011179238121*log(1 + exp(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3]))) + 0.7537599741964742*log(1 + exp(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))))) + 0.008349755673621395*log(1 + exp(-0.40141742757451926 + 0.7289092837988314*log(1 + exp(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3]))) + 0.16077592157977172*log(1 + exp(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3]))) + 0.6378484913808045*log(1 + exp(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3]))) + 0.9326240829024273*log(1 + exp(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3]))) + 0.9638370498324353*log(1 + exp(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))))) + -0.13663535205769417*log(1 + exp(0.17680227033398044 + 0.8293506390949057*log(1 + exp(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3]))) + -0.09356004658531791*log(1 + exp(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3]))) + -0.16358082594109513*log(1 + exp(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3]))) + 0.9029174038261001*log(1 + exp(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3]))) + 0.7318887196209252*log(1 + exp(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using Documenter, BibInternal

makedocs(
    sitename = "BibInternal.jl",
    authors = "Jean-François BAFFIER",
    repo="https://github.com/Humans-of-Julia/BibInternal.jl/blob/{commit}{path}#L{line}",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    pages = [
        "Entries" => "index.md",
        "BibTeX" => "bibtex.md",
        "Utilities" => "utilities.md"
    ]
)

deploydocs(; repo = "github.com/Humans-of-Julia/BibInternal.jl.git", devbranch = "master")

ext: jl


text:
module POMDPGifs

using Reel
using POMDPs
using POMDPSimulators
using POMDPModelTools
using POMDPPolicies
using ProgressMeter
using Parameters
using Random

export
    GifSimulator,
    makegif

struct SavedGif
    filename::String
end

function Base.show(io::IO, m::MIME"text/html", g::SavedGif)
    println(io, "<img src=\"$(g.filename)\">")
end

"""
    GifSimulator(<keyword arguments>)

Create a simulator for producing a gif output by calling `POMDPModelTools.render` at each step.

# Keyword Arguments
- `filename::String=tempname()*".gif"`
- `fps::Int=2`: frames per second
- `spec::Any`: specification for which elements of a step to render (see `POMDPSimulators.eachstep`)
- `max_steps::Int=nothing`
- `rng::AbstractRNG=GLOBAL_RNG`
- `show_progress::Bool`
- `extra_initial::Bool` if set to true, the simulator adds an extra step at time 0 (before first transition)
- `extra_final::Boll` if set to true, the simulator adds an extra setp at the end (after the last transition)
- `render_kwargs`: keyword args to be fed to `POMDPModelTools.render`
"""
@with_kw mutable struct GifSimulator <: Simulator
    filename::String                = tempname()*".gif"
    fps::Int                        = 2
    spec::Union{Nothing, Any}       = nothing
    max_steps::Union{Nothing, Int}  = nothing
    rng::AbstractRNG                = Random.GLOBAL_RNG
    show_progress::Bool             = max_steps != nothing
    extra_initial::Bool             = false
    extra_final::Bool               = false
    render_kwargs                   = NamedTuple()
end

function POMDPs.simulate(s::GifSimulator, m::Union{MDP, POMDP}, p::Policy=RandomPolicy(m, rng=s.rng), args...)

    # run simulation
    sim = HistoryRecorder(rng = s.rng,
                          max_steps = s.max_steps,
                          show_progress = s.show_progress
                         )
    hist = simulate(sim, m, p, args...)

    makegif(m, hist,
            filename=s.filename,
            spec=s.spec,
            show_progress=s.show_progress,
            extra_initial=s.extra_initial,
            extra_final=s.extra_final,
            render_kwargs=s.render_kwargs,
            fps=s.fps
           )
end

"""
    makegif(m; kwargs...)
    makegif(m, policy; kwargs...)
    makegif(m, policy, args...; kwargs...)

Create a gif of a single simulation of a POMDP or MDP by calling `POMDPModelTools.render` at each step.

# Arguments
- `m::Union{POMDP,MDP}`: the model to be simulated

All other positional arguments, for instance a policy, updater, initial state, etc. will be fed to the `simulate` function. See that documentation for more info.

# Keyword Arguments
All keyword arguments are fed to the `GifSimulator` constructor. See its documentation for more info.
"""
function makegif(args...; kwargs...)
    sim = GifSimulator(;kwargs...)
    return simulate(sim, args...)
end

"""
    makegif(m, history; kwargs...)

Create a gif from a POMDP or MDP and a history by calling `POMDPModelTools.render` at each step.

# Arguments
- `m::Union{POMDP,MDP}`: domain model
- `history::POMDPSimulators.SimHistory`: history of states, actions, etc. for the gif

# Keyword Arguments
- `filename::String=tempname()*".gif"`
- `fps::Int=2`: frames per second
- `spec::Any`: specification for which elements of a step to render (see `POMDPSimulators.eachstep`)
- `show_progress::Bool`
- `extra_initial::Bool` if set to true, the simulator adds an extra step at time 0 (before first transition)
- `extra_final::Bool` if set to true, the simulator adds an extra setp at the end (after the last transition)
- `render_kwargs`: keyword args to be fed to `POMDPModelTools.render`
"""
function makegif(m::Union{POMDP, MDP}, hist::POMDPSimulators.SimHistory;
                 filename=tempname()*".gif",
                 spec=nothing,
                 show_progress::Bool=true,
                 extra_initial::Bool=false,
                 extra_final::Bool=false,
                 render_kwargs=NamedTuple(),
                 fps::Int=2
                )

    # deal with the spec
    if spec == nothing
        steps = eachstep(hist)
    else
        steps = eachstep(hist, spec)
    end

    if extra_initial
        first_step = first(steps)
        extra_init_step = (t=0, sp=get(first_step, :s, missing), bp=get(first_step, :b, missing))
        steps = vcat(extra_init_step, collect(steps))
    end
    if extra_final
        last_step = last(steps)
        extra_final_step = (t=length(steps)+1, s=get(last_step, :sp, missing), b=get(last_step, :bp, missing), done=true)
        steps = vcat(collect(steps), extra_final_step)
    end

    # create gif
    frames = Frames(MIME("image/png"), fps=fps)
    if show_progress
        p = Progress(length(steps), 0.1, "Rendering $(length(steps)) steps...")
    end

    for step in steps
        push!(frames, render(m, step; pairs(render_kwargs)...))
        if show_progress
            next!(p)
        end
    end
    if show_progress
        @info "Creating Gif..."
    end
    write(filename, frames)
    if show_progress
        @info "Done Creating Gif."
    end
    return SavedGif(filename)
end

end # module


text:
push!(LOAD_PATH,"../src/")
import Documenter, AircraftSystems
Documenter.makedocs(sitename="My Documentation"#,
#	modules = [ModulesHere]
)

max_stars_repo_name: byuflowlab/AircraftSystems

max_issues_repo_name: byuflowlab/AircraftSystems

max_forks_repo_name: byuflowlab/AircraftSystems


text:
function check_node(nd::EAGO.Script.NodeInfo, type::JuMP._Derivatives.NodeType , indx::Int, child::Vector{Int})
    flag = true
    if nd.nodetype != type
        flag = false
        return flag
    end
    if nd.index != indx
        flag = false
        return flag
    end
    for i in 1:length(child)
        if nd.children[i] != child[i]
            flag = false
            return flag
        end
    end
    return flag
end

@testset "Test Trace Script" begin
    function f1(x)
        return sin(3.0*x[1]) + x[2]
    end
    tape1 = EAGO.Script.trace_script(f1,2)
    @test check_node(tape1.nd[1], JuMP._Derivatives.VARIABLE, 1, [-1])
    @test check_node(tape1.nd[2], JuMP._Derivatives.VARIABLE, 2, [-1])
    @test check_node(tape1.nd[3], JuMP._Derivatives.VALUE, 1, [-2])
    @test check_node(tape1.nd[4], JuMP._Derivatives.CALL, 3, [3,1])
    @test check_node(tape1.nd[5], JuMP._Derivatives.CALLUNIVAR, 15, [4])
    @test check_node(tape1.nd[6], JuMP._Derivatives.CALL, 1, [5,2])
    @test tape1.const_values[1] == 3.0
    @test tape1.num_valued[1] == false
    @test tape1.num_valued[2] == false
    @test tape1.num_valued[3] == true
    @test tape1.num_valued[4] == false
    @test tape1.num_valued[5] == false
    @test tape1.num_valued[6] == false
    @test tape1.set_trace_count == 6
    @test tape1.const_count == 1

    function f2(x)
        z = abs(x[1])
        y = 1.3
        for i in 1:2
            y += i*z
        end
        return sin(3.0*x[2]) + y
    end
    tape2 = EAGO.Script.trace_script(f2,2)
    @test check_node(tape2.nd[1], JuMP._Derivatives.VARIABLE, 1, [-1])
    @test check_node(tape2.nd[2], JuMP._Derivatives.VARIABLE, 2, [-1])
    @test check_node(tape2.nd[3], JuMP._Derivatives.CALLUNIVAR, 3, [1])
    @test check_node(tape2.nd[4], JuMP._Derivatives.VALUE, 1, [-2])
    @test check_node(tape2.nd[5], JuMP._Derivatives.CALL, 3, [4,3])
    @test check_node(tape2.nd[6], JuMP._Derivatives.VALUE, 2, [-2])
    @test check_node(tape2.nd[7], JuMP._Derivatives.CALL, 1, [6,5])
    @test check_node(tape2.nd[8], JuMP._Derivatives.VALUE, 3, [-2])
    @test check_node(tape2.nd[9], JuMP._Derivatives.CALL, 3, [8,3])
    @test check_node(tape2.nd[10], JuMP._Derivatives.CALL, 1, [7,9])
    @test check_node(tape2.nd[11], JuMP._Derivatives.VALUE, 4, [-2])
    @test check_node(tape2.nd[12], JuMP._Derivatives.CALL, 3, [11,2])
    @test check_node(tape2.nd[13], JuMP._Derivatives.CALLUNIVAR, 15, [12])
    @test check_node(tape2.nd[14], JuMP._Derivatives.CALL, 1, [13,10])
    @test tape2.const_values[1] == 1.0
    @test tape2.const_values[2] == 1.3
    @test tape2.const_values[3] == 2.0
    @test tape2.const_values[4] == 3.0
    @test tape2.num_valued[1] == false
    @test tape2.num_valued[2] == false
    @test tape2.num_valued[3] == false
    @test tape2.num_valued[4] == true
    @test tape2.num_valued[5] == false
    @test tape2.num_valued[6] == true
    @test tape2.num_valued[7] == false
    @test tape2.num_valued[8] == true
    @test tape2.num_valued[9] == false
    @test tape2.num_valued[10] == false
    @test tape2.num_valued[11] == true
    @test tape2.num_valued[12] == false
    @test tape2.num_valued[13] == false
    @test tape2.num_valued[14] == false
    @test tape2.set_trace_count == 14
    @test tape2.const_count == 4

    #=
    nds = EAGO.Script.tape_to_list(tape2)
    @test nds[1].index == 1
    @test nds[1].children[1] == -1
    @test nds[13].index == 15
    @test nds[13].children[1] == 12
    =#

    function f3(x)
        z = abs(x[1])::Float64
        return z
    end
    tape3 = EAGO.Script.trace_script(f3,2)
    #=
    @test tape3.const_count == 0
    @test tape3.set_trace_count == 3
    @test check_node(tape3.nd[1], JuMP._Derivatives.VARIABLE, 1, [-1])
    @test check_node(tape3.nd[2], JuMP._Derivatives.VARIABLE, 2, [-1])
    @test check_node(tape3.nd[3], JuMP._Derivatives.CALLUNIVAR, 3, [1])
    =#
end


text:
# make your scripts automatically re-activate your project
cd(@__DIR__)
using Pkg; Pkg.activate("../../..")

# load packages
using ml_for_pt
using DelimitedFiles
using LaTeXStrings
using Plots
using Random
ENV["GKSwstype"]="nul"

# set path to save folder
save_folder = "./results/"

# set critical value of tuning parameter
p_crit = 1.0f0

# define parameter ranges
p_min = 0.1f0
p_max = 3.0f0
dp = 0.05f0
p_range = collect(p_min:dp:p_max)
p_range_LBC = collect(p_min-dp/2:dp:p_max+dp/2)

# construct probability distribution
distr,samples = MLP.jump_distr(p_crit=p_crit)

# supervised learning using analytical expression
# returns optimal predictions, indicator, and loss
pred_opt_SL, indicator_opt_SL, loss_opt_SL = MLP.get_indicators_SL_analytical(samples, distr, p_range, dp, p_min, p_max)

# prediction-based method using analytical expression
# returns optimal predictions, indicator, and loss
pred_opt_PBM, indicator_opt_PBM, loss_opt_PBM = MLP.get_indicators_PBM_analytical(samples, distr, p_range, dp)

# learning by confusion using analytical expression
# returns optimal indicator and loss
indicator_opt_LBC, loss_opt_LBC = MLP.get_indicators_LBC_analytical(samples, distr, p_range, p_range_LBC)

# save output in save_folder
open(save_folder*"p_range_SL.txt", "w") do io
  writedlm(io, p_range)
end

open(save_folder*"p_range_PBM.txt", "w") do io
  writedlm(io, p_range)
end

open(save_folder*"p_range_LBC.txt", "w") do io
  writedlm(io, p_range_LBC)
end

open(save_folder*"pred_opt_SL.txt", "w") do io
  writedlm(io, pred_opt_SL)
end

open(save_folder*"indicator_opt_SL.txt", "w") do io
  writedlm(io, indicator_opt_SL)
end

open(save_folder*"loss_opt_SL.txt", "w") do io
  writedlm(io, [loss_opt_SL])
end

open(save_folder*"pred_opt_PBM.txt", "w") do io
  writedlm(io, pred_opt_PBM)
end

open(save_folder*"indicator_opt_PBM.txt", "w") do io
  writedlm(io, indicator_opt_PBM)
end

open(save_folder*"loss_opt_PBM.txt", "w") do io
  writedlm(io, [loss_opt_PBM])
end

open(save_folder*"indicator_opt_LBC.txt", "w") do io
  writedlm(io, indicator_opt_LBC)
end

open(save_folder*"loss_opt_LBC.txt", "w") do io
  writedlm(io, loss_opt_LBC)
end

# plot and save results in save_folder
plt = plot(p_range, pred_opt_SL, dpi=300, c="black", label=false)
vline!([p_crit], label=L"p_{c}", c="red")
xlabel!(L"$p$")
ylabel!(L"$\hat{y}^{\mathrm{opt}}_{\mathrm{SL}}$")
savefig(save_folder*"pred_opt_SL.png")

plt = plot(p_range[2:end-1], indicator_opt_SL, dpi=300, c="black", label=false)
vline!([p_crit], label=L"p_{c}", c="red")
xlabel!(L"$p$")
ylabel!(L"$I^{\mathrm{opt}}_{\mathrm{SL}}$")
savefig(save_folder*"indicator_opt_SL.png")

plt = plot(p_range, pred_opt_PBM, dpi=300, c="black", label=false)
vline!([p_crit], label=L"p_{c}", c="red")
xlabel!(L"$p$")
ylabel!(L"$\hat{y}^{\mathrm{opt}}_{\mathrm{PBM}}$")
savefig(save_folder*"pred_opt_PBM.png")

plt = plot(p_range[2:end-1], indicator_opt_PBM, dpi=300, c="black", label=false)
vline!([p_crit], label=L"p_{c}", c="red")
xlabel!(L"$p$")
ylabel!(L"$I^{\mathrm{opt}}_{\mathrm{PBM}}$")
savefig(save_folder*"indicator_opt_PBM.png")

plt = plot(p_range_LBC, indicator_opt_LBC, dpi=300, c="black", label=false)
vline!([p_crit], label=L"p_{c}", c="red")
xlabel!(L"$p$")
ylabel!(L"$I^{\mathrm{opt}}_{\mathrm{LBC}}$")
savefig(save_folder*"indicator_opt_LBC.png")

plt = plot(p_range_LBC, loss_opt_LBC, dpi=300, c="black", label=false)
vline!([p_crit], label=L"p_{c}", c="red")
xlabel!(L"$p$")
ylabel!(L"$\mathcal{L}^{\mathrm{opt}}_{\mathrm{LBC}}$")
savefig(save_folder*"loss_opt_LBC.png")


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

using Markdown
using InteractiveUtils

# ╔═╡ 4b149ad2-5bca-11eb-0104-0be8dd786177
begin
	import Pkg
	Pkg.activate(Base.current_project())
	using ReinforcementLearning
	using Flux
	using Statistics
	using Plots
end

# ╔═╡ 48977b2e-5bc9-11eb-2ad7-958ca5fd2ecb
md"""
# Chapter 8.6 Trajectory Sampling

The general function `run(policy, env, stop_condition, hook)` is very flexible and powerful. However, we are not restricted to use it only. In this notebook, we'll see how to use part of the components provided in `ReinforcementLearning.jl` to finish some specific experiments.

First, let's define the environment mentioned in **Chapter 8.6**:
"""

# ╔═╡ 6bf4f104-5bca-11eb-2b4a-69c4ba7ffadc
begin
	mutable struct TestEnv <: AbstractEnv
		transitions::Array{Int, 3}
		rewards::Array{Float64, 3}
		reward_table::Array{Float64, 2}
		terminate_prob::Float64
		# cache
		s_init::Int
		s::Int
		reward::Float64
		is_terminated::Bool
	end

	function TestEnv(;ns=1000, na=2, b=1, terminate_prob=0.1,init_state=1)
		transitions = rand(1:ns, b, na, ns)
		rewards = randn(b, na, ns)
		reward_table = randn(na, ns)
		TestEnv(
			transitions,
			rewards,
			reward_table,
			terminate_prob,
			init_state,
			init_state,
			0.,
			false
		)
	end

	function (env::TestEnv)(a::Int)
		t = rand() < 0.1
		bᵢ = rand(axes(env.transitions, 1))

		env.is_terminated = t
		if t
			env.reward = env.reward_table[a, env.s]
		else
			env.reward = env.rewards[bᵢ, a, env.s]
		end

		env.s = env.transitions[bᵢ, a, env.s]

	end

	RLBase.state_space(env::TestEnv) = Base.OneTo(1:size(env.rewards, 3))
	RLBase.action_space(env::TestEnv) = Base.OneTo(1:size(env.rewards, 2))

	function RLBase.reset!(env::TestEnv)
		env.s = env.s_init
		env.reward = 0.0
		env.is_terminated = false
	end

	RLBase.is_terminated(env::TestEnv) = env.is_terminated
	RLBase.state(env::TestEnv) = env.s
	RLBase.reward(env::TestEnv) = env.reward
end

# ╔═╡ 989f834a-5bca-11eb-363a-5f4afce32520
md"""
Note that this environment is not described very clearly on the book. Part of the information are inferred from the [lisp source code](http://incompleteideas.net/book/code/sampling2.lisp).

!!! info
	Actually the lisp code is also not perfect, I spent a whole afternoon to figure out the code logic. So good luck if you also want to understand it.

The definitions above are just like any other environment we've defined before in previous chapters. Now we'll add an extra function to make it work for our **planning** purpose.
"""

# ╔═╡ 97be6cb8-5bcb-11eb-3dc7-1b8105b206d7
"""
Return all the possible next states and corresponding reward.
"""
function successors(env::TestEnv, s, a)
    S = @view env.transitions[:, a, s]
    R = @view env.rewards[:, a, s]
    zip(R, S)
end

# ╔═╡ 0c52737a-5bcd-11eb-1546-21885703f303
γ = 0.9

# ╔═╡ aa061636-5bce-11eb-1381-bb00fb913f9b
n_sweep=10

# ╔═╡ a466db90-5bcc-11eb-3385-e5117be33410
"""
Here we are only interested in the performance of Q
with env starting at state `1`. Note here we're calculating
the discounted reward.
"""
function eval_Q(Q, env;n_eval=100)
	R = 0.
	for _ in 1:n_eval
		reset!(env)
		i = 0
		while !is_terminated(env)
			a = Q(state(env)) |> argmax  # greedy
			env(a)
			R += reward(env) * γ^i
			i += 1
		end
	end
	R/n_eval
end

# ╔═╡ 9bf4aea2-5bcd-11eb-230f-cf63a009a0f9
"""
Calculate the expected gain.
"""
function gain(Q,env,s,a)
    p = env.terminate_prob
    r = env.reward_table[a, s]
    p * r + (1-p) * mean(r̄ + γ * maximum(Q(s′)) for (r̄, s′) in  successors(env, s, a))
end

# ╔═╡ b92ba1e2-5bcb-11eb-16cf-f15929c2697b
function sweep(;b = 1, ns=1000)

    na = 2

    α=1.0
    p = 0.1

    env= TestEnv(;ns=ns, na=na, b=b, terminate_prob=p)
    Q = TabularQApproximator(;n_state=ns, n_action=na, opt=Descent(α))

    i = 1
    vals = [eval_Q(Q, env)]
    for _ in 1:n_sweep
        for s in 1:ns
            for a in 1:na
                G = gain(Q,env,s,a)
                update!(Q, (s,a) =>  Q(s, a) - G)
                if i % 100 == 0
                    push!(vals, eval_Q(Q, env))
                end
                i += 1
            end
        end
    end
    vals
end

# ╔═╡ 86538552-5bce-11eb-2cbf-ff6e175ef421
function on_policy(;b = 1, ns=1000)

    na = 2

    α=1.0
    p = 0.1

    env= TestEnv(;ns=ns, na=na, b=b, terminate_prob=p)
    Q = TabularQApproximator(;n_state=ns, n_action=na, opt=Descent(α))

    i = 1
    vals = [eval_Q(Q, env)]

	explorer = EpsilonGreedyExplorer(0.1)
	for i in 1:(n_sweep * ns * na)
		is_terminated(env) && reset!(env)
		s = state(env)
		a = Q(s) |> explorer
		env(a)
		G = gain(Q, env, s, a)
		update!(Q, (s,a) => Q(s,a) - G)
		if i % 100 == 0
			push!(vals, eval_Q(Q, env))
		end
	end

    vals
end

# ╔═╡ f29c9f18-5bcd-11eb-05b5-c72249e02975
begin
	fig_8_8 = plot(legend=:bottomright, xlabel="Computation time", ylabel="Value of start state under greedy policy")
	for b in [1, 3, 10]
		plot!(fig_8_8, mean(sweep(;b=b) for _ in 1:200), label="uniform b=$b")
		plot!(fig_8_8, mean(on_policy(;b=b) for _ in 1:200), label="on-policy b=$b")
	end
	fig_8_8
end

# ╔═╡ 8cba1ec8-5bcf-11eb-32c3-f38dee82139a
begin
	fig_8_8_2 = plot(legend=:bottomright, xlabel="Computation time", ylabel="Value of start state under greedy policy")
	plot!(fig_8_8_2, mean(sweep(;ns=10_000) for _ in 1:200), label="uniform")
	plot!(fig_8_8_2, mean(on_policy(;ns=10_000) for _ in 1:200), label="on-policy")
	fig_8_8_2
end

# ╔═╡ Cell order:
# ╟─48977b2e-5bc9-11eb-2ad7-958ca5fd2ecb
# ╠═4b149ad2-5bca-11eb-0104-0be8dd786177
# ╠═6bf4f104-5bca-11eb-2b4a-69c4ba7ffadc
# ╟─989f834a-5bca-11eb-363a-5f4afce32520
# ╠═97be6cb8-5bcb-11eb-3dc7-1b8105b206d7
# ╠═0c52737a-5bcd-11eb-1546-21885703f303
# ╠═aa061636-5bce-11eb-1381-bb00fb913f9b
# ╠═a466db90-5bcc-11eb-3385-e5117be33410
# ╠═9bf4aea2-5bcd-11eb-230f-cf63a009a0f9
# ╠═b92ba1e2-5bcb-11eb-16cf-f15929c2697b
# ╠═86538552-5bce-11eb-2cbf-ff6e175ef421
# ╠═f29c9f18-5bcd-11eb-05b5-c72249e02975
# ╠═8cba1ec8-5bcf-11eb-32c3-f38dee82139a


text:
function test_vec_interface()
  # serial tests only
  if MPI.Comm_size(MPI.COMM_WORLD) != 1
    return nothing
  end
  @testset "----- Testing Vec interface -----" begin
    const m = 3
    vals = PetscScalar[1, 2, 3]
    vals2 = zeros(vals)  # used with get_values
    idx = PetscInt[1, 2, 3]
    x_j = zeros(PetscScalar, 3)
    set_values1!(x_j, idx, vals)

    assembly_begin(x_j)
    assembly_end(x_j)

    # test set_values!, both modes
    @test isapprox( norm(x_j - vals), 0.0)
    set_values1!(x_j, idx, vals)
    @test isapprox( norm(x_j - vals), 0.0)
    set_values1!(x_j, idx, vals, ADD_VALUES)
    @test isapprox( norm(x_j - 2*vals), 0.0)
    get_values1!(x_j, idx, vals2)
    @test isapprox( norm(vals2 - 2*vals), 0.0)

    x_p = PetscVec(3, PETSc2.VECMPI, MPI.COMM_WORLD)
    set_values1!(x_p, idx, vals)
    assembly_begin(x_p)
    assembly_end(x_p)
    get_values1!(x_p, idx, vals2)
    @test isapprox( norm(vals2 - vals), 0.0)
    set_values1!(x_p, idx, vals)
    assembly_begin(x_p)
    assembly_end(x_p)
    get_values1!(x_p, idx, vals2)
    @test isapprox( norm(vals2 - vals), 0.0)
    set_values1!(x_p, idx, vals, ADD_VALUES)
    assembly_begin(x_p)
    assembly_end(x_p)
    get_values1!(x_p, idx, vals2)
    @test isapprox( norm(vals2 - 2*vals), 0.0)

    PetscView(x_p)  # if this doesn't error, that's all we can (easily) test

    @test ( size_local(x_j) )== (3,)
    @test ( size_global(x_j) )== (3,)
    @test ( size_local(x_p) )== (3,)
    @test ( size_global(x_p) )== (3,)
    @test ( length_local(x_j) )== 3
    @test ( length_global(x_j) )== 3
    @test ( length_local(x_p) )== 3
    @test ( length_global(x_p) )== 3
    @test ( local_indices(x_j) )== 1:3
    @test ( local_indices(x_p) )== 1:3

    fill_zero!(x_j)
    @test isapprox( norm(x_j), 0.0)
    fill_zero!(x_p)
    @test isapprox( norm(x_p), 0.0)

    fill!(x_j, 1.0)
    @test isapprox( norm(x_j, Inf), 1.0)
    fill!(x_p, 1.0)
    @test isapprox( norm(x_p, Inf), 1.0)

    set_values1!(x_j, idx, vals)
    set_values1!(x_p, idx, vals)
    assembly_begin(x_p)
    assembly_end(x_p)

    @test isapprox( maximum(x_p), 3.0)
    @test isapprox( minimum(x_p), 1.0)

    x2_p = copy(x_p)
    get_values1!(x2_p, idx, vals2)
    @test isapprox( norm(vals2 - vals), 0.0)

    fill_zero!(x2_p)
    get_values1!(x2_p, idx, vals2)
    @test isapprox( norm(vals2), 0.0)

    copy!(x2_p, x_p)
    get_values1!(x2_p, idx, vals2)
    @test isapprox( norm(vals2 - vals), 0.0)

    # diagonal_shift
    x_j[:] = vals
    alpha = 2
    diagonal_shift!(x_j, alpha)
    @test isapprox( norm(x_j - (vals + alpha)), 0.0) atol=1e-14

    set_values1!(x_p, idx, vals, INSERT_VALUES)
    assembly_begin(x_p); assembly_end(x_p)
    diagonal_shift!(x_p, alpha)
    get_values1!(x_p, idx, vals2)
    @test isapprox( norm(vals2 - (vals + alpha)), 0.0) atol=1e-14

    PetscDestroy(x_p)
    PetscDestroy(x2_p)
  end
end

test_vec_interface()


text:
__precompile__()

module FixedPointNumbers

using Base: reducedim_initarray

import Base: ==, <, <=, -, +, *, /, ~, isapprox,
             convert, promote_rule, show, showcompact, isinteger, abs, decompose,
             isnan, isinf, isfinite,
             zero, oneunit, one, typemin, typemax, realmin, realmax, eps, sizeof, reinterpret,
             float, trunc, round, floor, ceil, bswap,
             div, fld, rem, mod, mod1, fld1, min, max, minmax,
             start, next, done, rand
if isdefined(Base, :rem1)
    import Base: rem1
end
using Base: @pure

# T => BaseType
# f => Number of Bytes reserved for fractional part
abstract type FixedPoint{T <: Integer, f} <: Real end

export
    FixedPoint,
    Fixed,
    Normed,
# "special" typealiases
    # Q and U typealiases are exported in separate source files
# literal constructor constants
    uf8,
    uf10,
    uf12,
    uf14,
    uf16,
# Functions
    scaledual

reinterpret(x::FixedPoint) = x.i
reinterpret(::Type{T}, x::FixedPoint{T,f}) where {T,f} = x.i

# construction using the (approximate) intended value, i.e., N0f8
*(x::Real, ::Type{X}) where {X<:FixedPoint} = X(x)

# comparison
==(x::T, y::T) where {T <: FixedPoint} = x.i == y.i
 <(x::T, y::T) where {T <: FixedPoint} = x.i  < y.i
<=(x::T, y::T) where {T <: FixedPoint} = x.i <= y.i
"""
    isapprox(x::FixedPoint, y::FixedPoint; rtol=0, atol=max(eps(x), eps(y)))

For FixedPoint numbers, the default criterion is that `x` and `y` differ by no more than `eps`, the separation between adjacent fixed-point numbers.
"""
function isapprox(x::T, y::T; rtol=0, atol=max(eps(x), eps(y))) where {T <: FixedPoint}
    maxdiff = T(atol+rtol*max(abs(x), abs(y)))
    rx, ry, rd = reinterpret(x), reinterpret(y), reinterpret(maxdiff)
    abs(signed(widen1(rx))-signed(widen1(ry))) <= rd
end
function isapprox(x::FixedPoint, y::FixedPoint; rtol=0, atol=max(eps(x), eps(y)))
    isapprox(promote(x, y)...; rtol=rtol, atol=atol)
end

# predicates
isinteger(x::FixedPoint{T,f}) where {T,f} = (x.i&(1<<f-1)) == 0

# traits
typemax(::Type{T}) where {T <: FixedPoint} = T(typemax(rawtype(T)), 0)
typemin(::Type{T}) where {T <: FixedPoint} = T(typemin(rawtype(T)), 0)
realmin(::Type{T}) where {T <: FixedPoint} = eps(T)
realmax(::Type{T}) where {T <: FixedPoint} = typemax(T)

widen1(::Type{Int8})   = Int16
widen1(::Type{UInt8})  = UInt16
widen1(::Type{Int16})  = Int32
widen1(::Type{UInt16}) = UInt32
widen1(::Type{Int32})  = Int64
widen1(::Type{UInt32}) = UInt64
widen1(::Type{Int64})  = Int128
widen1(::Type{UInt64}) = UInt128
widen1(::Type{UInt128}) = UInt128
widen1(x::Integer) = x % widen1(typeof(x))

const ShortInts = Union{Int8,UInt8,Int16,UInt16}

floattype(::Type{FixedPoint{T,f}}) where {T <: ShortInts,f} = Float32
floattype(::Type{FixedPoint{T,f}}) where {T <: Integer,f} = Float64
floattype(::Type{F}) where {F <: FixedPoint} = floattype(supertype(F))
floattype(x::FixedPoint) = floattype(typeof(x))

nbitsfrac(::Type{FixedPoint{T,f}}) where {T <: Integer,f} = f
nbitsfrac(::Type{F}) where {F <: FixedPoint} = nbitsfrac(supertype(F))

rawtype(::Type{FixedPoint{T,f}}) where {T <: Integer,f} = T
rawtype(::Type{F}) where {F <: FixedPoint} = rawtype(supertype(F))
rawtype(x::FixedPoint) = rawtype(typeof(x))

# This IOBuffer is used during module definition to generate typealias names
_iotypealias = IOBuffer()

# Printing. These are used to generate type-symbols, so we need them
# before we include any files.
function showtype(io::IO, ::Type{X}) where {X <: FixedPoint}
    print(io, typechar(X))
    f = nbitsfrac(X)
    m = sizeof(X)*8-f-signbits(X)
    print(io, m, 'f', f)
    io
end
function show(io::IO, x::FixedPoint{T,f}) where {T,f}
    showcompact(io, x)
    showtype(io, typeof(x))
end
const _log2_10 = 3.321928094887362
showcompact(io::IO, x::FixedPoint{T,f}) where {T,f} = show(io, round(convert(Float64,x), ceil(Int,f/_log2_10)))

if VERSION >= v"0.7.0-DEV.1790"
    function Base.showarg(io::IO, a::Array{T}, toplevel) where {T<:FixedPoint}
        toplevel || print(io, "::")
        print(io, "Array{")
        showtype(io, T)
        print(io, ",$(ndims(a))}")
        toplevel && print(io, " with eltype ", T)
    end
end

include("fixed.jl")
include("normed.jl")
include("deprecations.jl")

eps(::Type{T}) where {T <: FixedPoint} = T(oneunit(rawtype(T)),0)
eps(::T) where {T <: FixedPoint} = eps(T)
sizeof(::Type{T}) where {T <: FixedPoint} = sizeof(rawtype(T))

# Promotions for reductions
const Treduce = Float64
if isdefined(Base, :r_promote)
    # Julia v0.6
    Base.r_promote(::typeof(+), x::FixedPoint{T}) where {T} = Treduce(x)
    Base.r_promote(::typeof(*), x::FixedPoint{T}) where {T} = Treduce(x)
    Base.reducedim_init(f::typeof(identity),
                   op::typeof(+),
                   A::AbstractArray{T}, region) where {T <: FixedPoint} =
                       Base.reducedim_initarray(A, region, zero(Treduce))
    Base.reducedim_init(f::typeof(identity),
                   op::typeof(*),
                   A::AbstractArray{T}, region) where {T <: FixedPoint} =
                       Base.reducedim_initarray(A, region, oneunit(Treduce))
else
    # Julia v0.7
    Base.add_sum(x::FixedPoint, y::FixedPoint) = Treduce(x) + Treduce(y)
    Base.reduce_empty(::typeof(Base.add_sum), ::Type{F}) where {F<:FixedPoint}  = zero(Treduce)
    Base.reduce_first(::typeof(Base.add_sum), x::FixedPoint)   = Treduce(x)
    Base.mul_prod(x::FixedPoint, y::FixedPoint) = Treduce(x) * Treduce(y)
    Base.reduce_empty(::typeof(Base.mul_prod), ::Type{F}) where {F<:FixedPoint} = one(Treduce)
    Base.reduce_first(::typeof(Base.mul_prod), x::FixedPoint)  = Treduce(x)
end

for f in (:div, :fld, :fld1)
    @eval begin
        $f(x::T, y::T) where {T <: FixedPoint} = $f(reinterpret(x),reinterpret(y))
    end
end
for f in (:rem, :mod, :mod1, :rem1, :min, :max)
    if f === :rem1 && !isdefined(Base, :rem1)
        continue
    end
    @eval begin
        $f(x::T, y::T) where {T <: FixedPoint} = T($f(reinterpret(x),reinterpret(y)),0)
    end
end

# When multiplying by a float, reduce two multiplies to one.
# Particularly useful for arrays.
scaledual(Tdual::Type, x) = oneunit(Tdual), x
scaledual(b::Tdual, x) where {Tdual <: Number} = b, x
scaledual(Tdual::Type, x::Union{T,AbstractArray{T}}) where {T <: FixedPoint} =
    convert(Tdual, 1/oneunit(T)), reinterpret(rawtype(T), x)
scaledual(b::Tdual, x::Union{T,AbstractArray{T}}) where {Tdual <: Number,T <: FixedPoint} =
    convert(Tdual, b/oneunit(T)), reinterpret(rawtype(T), x)

@noinline function throw_converterror(::Type{T}, x) where {T <: FixedPoint}
    n = 2^(8*sizeof(T))
    bitstring = sizeof(T) == 1 ? "an 8-bit" : "a $(8*sizeof(T))-bit"
    io = IOBuffer()
    showcompact(io, typemin(T)); Tmin = String(take!(io))
    showcompact(io, typemax(T)); Tmax = String(take!(io))
    throw(ArgumentError("$T is $bitstring type representing $n values from $Tmin to $Tmax; cannot represent $x"))
end

rand(::Type{T}) where {T <: FixedPoint} = reinterpret(T, rand(rawtype(T)))
rand(::Type{T}, sz::Dims) where {T <: FixedPoint} = reinterpret(T, rand(rawtype(T), sz))

end # module


text:
"""
General observers
"""

abstract type AbstractObserver{T} end

"""
    LinearObserver(H::Matrix[, R::Matrix])
    LinearObserver(H::AbstractFloat[, R::AbstractFloat])
    LinearObserver(H::Vector[, R::Matrix])
    LinearObserver(H::RowVector[, R::Union{Matrix, AbstractFloat}])

Linear observer with the form

\$y_k = H x_k + v_k \\quad \\text{where} \\quad  v_k \\sim N(0, R)\$

A linear observer can be constructed with both matrix and scalar inputs. The
observer will not model measurement noise if the covariance R is not provided.
"""
struct LinearObserver{T<:AbstractFloat} <: AbstractObserver{T}
    H::Matrix{T}
    R::Covariance{T}

    function LinearObserver(H::Matrix{T}, R::Covariance{T}) where T
        if (size(H,1), size(H,1)) != size(R)
            throw(DimensionMismatch("Incompatible size of system matrices."))
        end
        new{T}(H, R)
    end
end
LinearObserver{T}(H::Matrix{T}) = LinearObserver(H,zeros(T,size(H,1),size(H,1)))
LinearObserver{T<:AbstractFloat}(H::T) =
    LinearObserver(reshape([H],1,1), zeros(T,1,1))
LinearObserver{T<:AbstractFloat}(H::T, R::T) =
    LinearObserver(reshape([H],1,1), reshape([R],1,1))
LinearObserver{T<:AbstractFloat}(H::Matrix{T}, R::T) =
    LinearObserver(H, reshape([R],1,1))

LinearObserver(H::Vector) = LinearObserver(reshape(H,length(H),1))
LinearObserver{T<:AbstractFloat}(H::Vector{T}, R::Covariance{T}) =
    LinearObserver(reshape(H,length(H),1), R)

LinearObserver(H::RowVector) = LinearObserver(Array(H))
LinearObserver{T<:AbstractFloat}(H::RowVector{T}, R::T) =
    LinearObserver(Array(H), reshape([R],1,1))
LinearObserver{T<:AbstractFloat}(H::RowVector{T}, R::Covariance{T}) =
    LinearObserver(Array(H), R)

"""
    NonlinearObserver(H::Function, dH_dx::Function, R::Matrix)
    NonlinearObserver(H::Function, dH_dx::Function, R::AbstractFloat)

Nonlinear observer with the form

\$y_k = H(k, x_k) + v \\quad  \\text{where} \\quad v \\sim N(0, R)\$

NonlinearObserver constructors require both the measurement function,
`H(t::Number, x::Vector)`, and its Jacobian, `dH_dx(t::Number, x::Vector)`. Both
`H()` and `dH_dx()` must return a vector.
"""
struct NonlinearObserver{T<:AbstractFloat} <: AbstractObserver{T}
    H::Function
    dH_dx::Function
    R::Covariance{T}

    function NonlinearObserver(H::Function, dH_dx::Function,
                               R::Covariance{T}) where T
        new{T}(H, dH_dx, R)
    end
end
NonlinearObserver(H::Function, dH_dx::Function, R::AbstractFloat) =
    NonlinearObserver(H, dH_dx, reshape([R], 1, 1))
function NonlinearObserver(lin_obs::LinearObserver)
    H(t,x) = lin_obs.H*x
    dH_dx(t,x) = lin_obs.H
    return NonlinearObserver(H, dH_dx, lin_obs.R)
end

"""
    assert_compatibility(obs::LinearObserver{T}, state::AbstractState{T})

Require linear observer dimensions to be compatible with input state.

    assert_compatibility(state::AbstractState{T}, obs::LinearObserver{T})

Require linear observer dimensions to be left-compatible with input state.
"""
function assert_compatibility{T}(obs::LinearObserver{T},state::AbstractState{T})
    if size(obs.H, 2) != length(state.x)
        throw(DimensionMismatch("Observer incompatible with input state."))
    end
end
function assert_compatibility{T}(state::AbstractState{T},obs::AbstractObserver{T})
    if size(obs.R, 1) != length(state.x)
        throw(DimensionMismatch("Observer incompatible with input state."))
    end
end

"""
    observable(A::Matrix, H::Matrix)
    observable(sys::LinearSystem, obs::LinearObserver)

Evaluate the observability of a linear model. Returns a boolean True/False
result.
"""
function observable(A, H)
    n = size(A, 2)
    p = size(H, 1)
    observability = zeros(n*p, n)

    observability[1:p, :] .= H
    for idx = 1:n-1
        observability[p*idx+1:p*(idx+1), :] .= H * A
    end

    rank(observability) == n ? true : false
end
observable(sys::LinearSystem, obs::LinearObserver) = observable(sys.A, obs.H)

"""
    predict(state::AbstractState, obs::AbstractObserver)

Predict a state through an arbitrary observer. The `predict` method will advance
both the state vector and, if necessary, its covariance through the input
measurement model. For NonlinearObservers, the covariance prediction is only a
linear approximation.
"""
function predict(::AbstractObserver) end
function predict{T}(state::DiscreteState{T}, obs::LinearObserver{T})
    return DiscreteState(obs.H*state.x, state.t)
end
function predict{T}(state::ContinuousState{T}, obs::LinearObserver{T})
    return ContinuousState(obs.H*state.x, state.t)
end
function predict{T}(state::UncertainDiscreteState{T}, obs::LinearObserver{T})
    return UncertainDiscreteState(
        obs.H*state.x, obs.H*state.P*obs.H'+obs.R, state.t)
end
function predict{T}(state::UncertainContinuousState{T}, obs::LinearObserver{T})
    return UncertainContinuousState(
        obs.H*state.x, obs.H*state.P*obs.H'+obs.R, state.t)
end
function predict{T}(state::DiscreteState{T}, obs::NonlinearObserver{T})
    return DiscreteState(obs.H(state.t, state.x), state.t)
end
function predict{T}(state::ContinuousState{T}, obs::NonlinearObserver{T})
    return ContinuousState(obs.H(state.t, state.x), state.t)
end
function predict{T}(state::UncertainDiscreteState{T}, obs::NonlinearObserver{T})
    jac = obs.dH_dx(state.t, state.x)
    return UncertainDiscreteState(
        obs.H(state.t, state.x), jac*state.P*jac'+obs.R, state.t)
end
function predict{T}(state::UncertainContinuousState{T},
                    obs::NonlinearObserver{T})
    jac = obs.dH_dx(state.t, state.x)
    return UncertainContinuousState(
        obs.H(state.t, state.x), jac*state.P*jac'+obs.R, state.t)
end

"""
   simulate(state::AbstractState, obs::AbstractObserver)

Simulate a state observation with initial state error and process noise. This
method returns an absolute state by sampling the initial state (if uncertain),
and propagating that state through the input system in the presence of
measurement noise.
"""
function simulate(::AbstractObserver) end
function simulate(state::AbstractAbsoluteState, obs::AbstractObserver)
    return simulate(make_uncertain(state), obs)
end
function simulate(state::AbstractUncertainState, obs::LinearObserver)
    return sample(predict(state, obs))
end
function simulate(state::AbstractUncertainState, obs::NonlinearObserver)
    return sample(predict(make_uncertain(sample(state)), obs))
end


text:
@testset "Conjugacy" begin
  include("Conjugacy/Conjugacy.jl")
  include("Conjugacy/Husert.jl")
end


text:
using GraphPlot

"""
    drawdag(dag; layout=spring_layout, ...)

Draws dag with given layout method (check GraphPlot package for their description)
"""
drawdag(dag::DAG; layout::Function=spring_layout, keyargs...) =
    drawdag(dag, layout(dag.graph)...; keyargs...)

"""
    drawdag(dag, locs_x, locs_y)

Draws dag with given locations of every node

# Examples

```julia
julia> g = DAG(:A => :C, :C => :B)
DAG: {3, 2} directed simple Int64 graph with labels [:A, :B, :C])

julia> drawdag(g, [0, 0, 1], [0, 1, 1])
```
"""
function drawdag(dag::DAG, locs_x, locs_y; keyargs...)
    gplot(dag.graph, locs_x, locs_y; nodelabel=dag.labels)
end

export
    drawdag


text:
using Clingo
x,y = Clingo.testsolve("#const n=4. c(1..n). 1 {color(X,I) : c(I)} 1 :- v(X). :- color(X,I), color(Y,I), e(X,Y), c(I). v(1..100). e(1,55).")


text:
#
# Hierarchical closure.
#

function ConnectQuery(refl::Bool, q::Query)
    @assert fits(output(q), input(q))
    @assert isoptional(output(q))
    Query(
        ConnectSig(refl),
        [q],
        input(q),
        output(q) |> setoptional(!refl) |> setplural())
end

immutable ConnectSig <: AbstractSignature
    refl::Bool
end

ev(sig::ConnectSig, args::Vector{Query}, ity::Input, oty::Output, iflow::InputFlow) =
    ev(sig, args..., ity, oty, iflow)

function ev(sig::ConnectSig, arg::Query, ::Input, oty::Output, iflow::InputFlow)
    cols = Column[]
    if sig.refl
        push!(cols, Column(OneTo(length(iflow)+1), values(iflow)))
    end
    oflow = ev(arg, iflow)
    while length(values(oflow)) > 0
        push!(cols, column(oflow))
        iflow = distribute(iflow, oflow)
        oflow = ev(arg, iflow)
    end
    if isempty(cols)
        pile = column(oflow)
    else
        pile = cols[end]
        for k = endof(cols)-1:-1:1
            pile = connect_impl(cols[k], pile)
        end
    end
    return OutputFlow(oty, pile)
end

function connect_impl(col1::Column, col2::Column)
    len1 = length(col1)
    len2 = length(col2)
    offs1 = offsets(col1)
    offs2 = offsets(col2)
    vals1 = values(col1)
    vals2 = values(col2)
    L = length(vals1)
    offs = Vector{Int}(len1+1)
    offs[1] = 1
    vals = vcat(vals1, vals2)
    idxs = Vector{Int}(length(vals))
    n = 1
    for i = 1:len1
        l1 = offs1[i]
        r1 = offs1[i+1]
        for j = l1:r1-1
            idxs[n] = j
            n += 1
            l2 = offs2[j]
            r2 = offs2[j+1]
            for k = l2:r2-1
                idxs[n] = L + k
                n += 1
            end
        end
        offs[i+1] = n
    end
    return Column(offs, vals[idxs])
end

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
module backend

using PyCall
include("utils.jl")

export set_backend,
    context,
    tensor,
    is_tensor,
    shape,
    ndim,
    to_numpy,
    copy,
    transpose,
    ones,
    ones_like,
    zeros,
    zeros_like,
    sum,
    norm,
    dot,
    power,
    array_equal,
    einsum,
    random,
    seed,
    tensorinv

const T = PyNULL()

function __init__()
    copy!(T, pyimport_conda("autohoot.backend", "autohoot"))
end

@func_loader T set_backend context tensor is_tensor shape ndim to_numpy copy transpose ones ones_like zeros zeros_like
@func_loader T sum norm dot power array_equal einsum random seed tensorinv

end


text:
module config

import Hwloc

"Get the number of available CPUs."
function available_cpu_count()
    topology = Hwloc.topology_load();
    counts = Hwloc.histmap(topology);
    return counts[:Core];
end

"Get the number of available PUs."
function available_pu_count()
    topology = Hwloc.topology_load();
    counts = Hwloc.histmap(topology);
    return counts[:PU];
end

end  # module


text:
# aux functions
ns(l) = Int(net[:n_s][l])
nr(l) = Int(net[:n_r][l])
z(x) = quantile(Normal(0,1),1-x)

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table! s begin
        "--α_ϱ", "-p"
            help = "pressure penalty α_ϱ"
            arg_type = Float64
            default = 0.0
        "--α_ϑ", "-i"
            help = "injection variability coefficient α_ϑ"
            arg_type = Float64
            default = 0.025
        "--α_ψ", "-l"
            help = "linepack variability coefficient α_ψ"
            arg_type = Float64
            default = 100.0
        "--σ²", "-s"
            help = "wind covariance"
            arg_type = Float64
            default = 0.15
        "--ε̅_p", "-e"
            help = "electric power constraint violation tolerance"
            arg_type = Float64
            default = 0.01
        "--ε̅_g", "-g"
            help = "gas constraint violation tolerance"
            arg_type = Float64
            default = 0.005
        "--Ns", "-n"
            help = "number of samples for out-of-sample analysis"
            arg_type = Int64
            default = 1000
    end
    return parse_args(s)
end
args = parse_commandline()
outdir = "output"
mkpath(outdir)

function remove_col_and_row(B,refbus)
    @assert size(B,1) == size(B,2)
    n = size(B,1)
    return B[1:n .!= refbus, 1:n .!= refbus]
end

function full_matrix(A,ref)
    Nb = size(A,1)+1
    V = zeros(Nb,Nb)
    for i in 1:Nb, j in 1:Nb
        i < ref && j < ref ? V[i,j] = A[i,j] : NaN
        i > ref && j > ref ? V[i,j] = A[i-1,j-1] : NaN
        i > ref && j < ref ? V[i,j] = A[i-1,j] : NaN
        i < ref && j > ref ? V[i,j] = A[i,j-1] : NaN
    end
    return V
end

function build_B̆(B̂inv,refbus)
    Nb = size(B̂inv,1)+1
    B̆ = zeros(Nb,Nb)
    for i in 1:Nb, j in 1:Nb
        if i < refbus && j < refbus
            B̆[i,j] = B̂inv[i,j]
        end
        if i > refbus && j > refbus
            B̆[i,j] = B̂inv[i-1,j-1]
        end
        if i > refbus && j < refbus
            B̆[i,j] = B̂inv[i-1,j]
        end
        if i < refbus && j > refbus
            B̆[i,j] = B̂inv[i,j-1]
        end
    end
    return B̆
end

function W(i)
    W = zeros(length(ξ[:bus]),set[:k_t][i])
    for t in 1:i
        W[:,1] .= ξ[:w_f][1] * set[:w̅]
        if t > 1
            for i in 1:length(ξ[:bus])
                if i <= 4
                    W[i,set[:k_t][t]-2] = (ξ[:w_f][t] - ξ[:w_f][t-1]) * set[:w̅]
                end
                if i > 4 && i <= 9
                    W[i,set[:k_t][t]-1] = (ξ[:w_f][t] - ξ[:w_f][t-1]) * set[:w̅]
                end
                if i > 9
                    W[i,set[:k_t][t]] = (ξ[:w_f][t] - ξ[:w_f][t-1]) * set[:w̅]
                end
            end
        end
    end
    return W
end

function D(i)
    D = zeros(net[:D],set[:k_t][i])
    for t in 1:i
        D[:,1] .= net[:d]
    end
    return D
end

function S(t)
    S = zeros(set[:k_t][t],set[:k_t][set[:T]])
    S[diagind(S)] .= 1
    return S
end

function Δ(i)
    L = zeros(length(net[:N]),set[:k_t][i])
    L[:,1] = net[:δ]

    L_opf = CSV.File("data/gas_uncertainty.csv",header=0) |> Tables.matrix
    # L_opf = [zeros(i)';L_opf].*10
    L_opf = L_opf[:,1:set[:k_t][i]] * 3
    L[findall(x->x>0,L_opf[:,1]),:] = L_opf[findall(x->x>0,L_opf[:,1]),:]
    return L
end


text:
# functions related to binomial distribution

import .RFunctions:
    binompdf,
    binomlogpdf,
    binomcdf,
    binomccdf,
    binomlogcdf,
    binomlogccdf,
    binominvcdf,
    binominvccdf,
    binominvlogcdf,
    binominvlogccdf

# pdf for numbers with generic types
binompdf(n::Real, p::Real, k::Real) = exp(binomlogpdf(n, p, k))

# logpdf for numbers with generic types
binomlogpdf(n::Real, p::Real, k::Real) = -log1p(n) - lbeta(n - k + 1, k + 1) + k * log(p) + (n - k) * log1p(-p)


text:
module Configuration

function notebook_path_suggestion()
    preferred_dir = startswith(Sys.BINDIR, pwd()) ? homedir() : pwd()
    return joinpath(preferred_dir, "") # so that it ends with / or \
end

"""
The HTTP server options. See `SecurityOptions` for additional settings.
"""
Base.@kwdef mutable struct ServerOptions
    root_url::Union{Nothing,String} = nothing
    host::String = "127.0.0.1"
    port::Union{Nothing,Integer} = nothing
    launch_browser::Bool = true
    show_file_system::Bool = true
    notebook_path_suggestion::String = notebook_path_suggestion()
    notebook::Union{Nothing,String} = nothing
end

"""
    SecurityOptions([; kwargs...])

Security settings for the HTTP server. Options are:

- `require_secret_for_open_links::Bool = true`

    Whether the links `http://localhost:1234/open?path=/a/b/c.jl`  and `http://localhost:1234/open?path=http://www.a.b/c.jl` should be protected.

    Use `true` for almost every setup. Only use `false` if Pluto is running in a safe container (like mybinder.org), where arbitrary code execution is not a problem.

- `require_secret_for_access::Bool = true`

    If false, you do not need to use a `secret` in the URL to access Pluto: you will be authenticated by visiting `http://localhost:1234/` in your browser. An authentication cookie is still used for access (to prevent XSS and deceptive links or an img src to `http://localhost:1234/open?url=badpeople.org/script.jl`), and is set automatically, but this request to `/` is protected by cross-origin policy.

    Use `true` on a computer used by multiple people simultaneously. Only use `false` if necessary.

**Leave these options on `true` for the most secure setup.**

Note that Pluto is quickly evolving software, maintained by designers, educators and enthusiasts — not security experts. If security is a serious concern for your application, then we recommend running Pluto inside a container and verifying the relevant security aspects of Pluto yourself.
"""
Base.@kwdef mutable struct SecurityOptions
    require_secret_for_open_links::Bool = true
    require_secret_for_access::Bool = true
end

"""
For internal use only.
"""
Base.@kwdef mutable struct EvaluationOptions
    run_notebook_on_load::Bool = true
    workspace_use_distributed::Bool = true
end

"""
These options will be passed as command line argument to newly launched processes.

The ServerSession contains a global version of this configuration, and each notebook can also have its own version.
"""
Base.@kwdef mutable struct CompilerOptions
    compile::Union{Nothing,String} = nothing
    sysimage::Union{Nothing,String} = nothing
    banner::Union{Nothing,String} = nothing
    optimize::Union{Nothing,Int} = nothing
    math_mode::Union{Nothing,String} = nothing

    # notebook specified options
    # the followings are different from
    # the default julia compiler options

    # we use nothing to represent "@v#.#"
    project::Union{Nothing,String} = "@."
    # we don't load startup file in notebook
    startup_file::Union{Nothing,String} = "no"
    # we don't load history file in notebook
    history_file::Union{Nothing,String} = "no"

    @static if VERSION > v"1.5.0-"
        threads::Union{Nothing,String} = nothing
    end
end # struct CompilerOptions

"""
Collection of all settings that configure a Pluto session.

`ServerSession` contains a `Configuration`.
"""
Base.@kwdef struct Options
    server::ServerOptions = ServerOptions()
    security::SecurityOptions = SecurityOptions()
    evaluation::EvaluationOptions = EvaluationOptions()
    compiler::CompilerOptions = CompilerOptions()
end

# We don't us an abstract type because Base.@kwdef does not support subtyping in Julia 1.0, only in ≥1.1
AbstractOptions = Union{EvaluationOptions,CompilerOptions,ServerOptions,SecurityOptions,Options}

function overlayed(original::AbstractOptions; changes...)
    new_kwargs = Dict()
    for name in fieldnames(typeof(original))
        new_kwargs[name] = get(changes, name, getfield(original, name))
    end
    return typeof(original)(;new_kwargs...)
end

# NOTE: printings are copy-pastable
function Base.show(io::IO, x::AbstractOptions)
    indent = get(io, :indent, 0)

    summary(io, x)
    println(io, "(")
    fnames = fieldnames(typeof(x))
    for each in fieldnames(typeof(x))
        print(IOContext(io, :indent => 2), " "^indent, " "^2, each, " = ", repr(getfield(x, each)))
        println(io, ", ")
    end
    print(io, " "^indent, ")")
    return
end

function from_flat_kwargs(; kwargs...)::Options
    server_options = Dict()
    security_options = Dict()
    evaluation_options = Dict()
    compiler_options = Dict()

    for (k, v) in kwargs
        if k in fieldnames(EvaluationOptions)
            evaluation_options[k] = v
        elseif k in fieldnames(CompilerOptions)
            compiler_options[k] = v
        elseif k in fieldnames(ServerOptions)
            server_options[k] = v
        elseif k in fieldnames(SecurityOptions)
            security_options[k] = v
        else
            throw(ArgumentError("""Key $k not recognised. Options are:\n$(join(
            [
                "Server Options:",
                map(x -> " "^2 * string(x), fieldnames(ServerOptions))...,
                "",
                "Security Options:",
                map(x -> " "^2 * string(x), fieldnames(SecurityOptions))...,
                "",
                "Evaluation Options:",
                map(x -> " "^2 * string(x), fieldnames(EvaluationOptions))...,
                "",
                "Compiler Options:",
                map(x -> " "^2 * string(x), fieldnames(CompilerOptions))...,
            ], '\n'))

            These can be used as keywords arguments to Pluto.run, or to Pluto.from_flat_kwargs to create a Pluto.Options object."""))
        end
    end

    return Options(
        server=ServerOptions(; server_options...),
        security=SecurityOptions(; security_options...),
        evaluation=EvaluationOptions(; evaluation_options...),
        compiler=CompilerOptions(; compiler_options...),
    )
end

end


text:
mutable struct ThermalGeneratorOperations{R,G,T,P}

    # Parameters

    fixedcost::Matrix{Float64}    # $/unit/investment period, g
    variablecost::Matrix{Float64} # $/MW/hour, r x g
    startupcost::Matrix{Float64}  # $/startup/unit, r x g
    shutdowncost::Matrix{Float64} # $/shutdown/unit, r x g

    # Variables

    committed::Array{VariableRef,4} # Units committed (r x g x t x p)
    started::Array{VariableRef,4}   # Units started up
    shutdown::Array{VariableRef,4}  # Units shut down

    energydispatch::Array{VariableRef,4} # Energy dispatch (MW, r x g x t x p)
    raisereserve::Array{VariableRef,4}   # Raise reserves provisioned
    lowerreserve::Array{VariableRef,4}   # Lower reserves provisioned

    # Expressions

    fixedcosts::Matrix{ExpressionRef} # ($, r x g)
    variablecosts::Array{ExpressionRef,3} # ($, r x g x p)
    operatingcosts::Matrix{ExpressionRef} # ($, r x g)

    ucap::Vector{ExpressionRef} # (MW, r)
    totalenergy::Array{ExpressionRef,3} # MW, r x t x p
    totalraisereserve::Array{ExpressionRef,3}
    totallowerreserve::Array{ExpressionRef,3}

    # Constraints

    minunitcommitment::Array{GreaterThanConstraintRef,4} # (r x g x t x p)
    maxunitcommitment::Array{LessThanConstraintRef,4}
    minunitstartups::Array{GreaterThanConstraintRef,4}
    maxunitstartups::Array{LessThanConstraintRef,4}   # (r x g x t x p)
    minunitshutdowns::Array{GreaterThanConstraintRef,4}  # (r x g x t x p)
    maxunitshutdowns::Array{LessThanConstraintRef,4}  # (r x g x t x p)

    unitcommitmentcontinuity::Array{EqualToConstraintRef,4}  # (r x g x t x p)
    minunituptime::Array{LessThanConstraintRef,4} # (r x g x t x p)
    minunitdowntime::Array{LessThanConstraintRef,4}

    mingeneration::Array{GreaterThanConstraintRef,4} # (r x g x t x p)
    maxgeneration::Array{LessThanConstraintRef,4}
    minlowerreserve::Array{GreaterThanConstraintRef,4}
    maxlowerreserve::Array{LessThanConstraintRef,4}
    minraisereserve::Array{GreaterThanConstraintRef,4}
    maxraisereserve::Array{LessThanConstraintRef,4}

    maxrampdown::Array{GreaterThanConstraintRef,4} # (r x g x t x p)
    maxrampup::Array{LessThanConstraintRef,4}

    function ThermalGeneratorOperations{T,P}(
        fixedcost::Matrix{Float64},
        variablecost::Matrix{Float64},
        startupcost::Matrix{Float64},
        shutdowncost::Matrix{Float64}
    ) where {T,P}

        R, G = size(fixedcost)
        @assert size(variablecost) == (R, G)
        @assert size(startupcost) == (R, G)
        @assert size(shutdowncost) == (R, G)

        new{R,G,T,P}(fixedcost, variablecost, startupcost, shutdowncost)

    end

end

function ThermalGeneratorOperations{T,P}(
    thermalgens::ThermalGenerators{G}, regionlookup::Dict{String,Int},
    thermalpath::String
) where {G,T,P}

    R = length(regionlookup)
    thermaldata = DataFrame!(CSV.File(joinpath(thermalpath, "parameters.csv"),
                                      types=scenarios_resource_param_types))

    thermallookup = Dict(zip(thermalgens.name, 1:G))
    fixedcost = zeros(Float64, R, G)
    variablecost = zeros(Float64, R, G)
    startupcost = zeros(Float64, R, G)
    shutdowncost = zeros(Float64, R, G)

    for row in eachrow(thermaldata)

        gen_idx = thermallookup[row.class]
        region_idx = regionlookup[row.region]

        fixedcost[region_idx, gen_idx] = row.fixedcost
        variablecost[region_idx, gen_idx] = row.variablecost
        startupcost[region_idx, gen_idx] = row.startupcost
        shutdowncost[region_idx, gen_idx] = row.shutdowncost

    end

    return ThermalGeneratorOperations{T,P}(
        fixedcost, variablecost, startupcost, shutdowncost)

end

function setup!(
    ops::ThermalGeneratorOperations{R,G,T,P},
    units::ThermalGenerators{G},
    m::Model, invs::ResourceInvestments{R,G},
    periodweights::Vector{Float64},
    s::AbstractScenario
) where {R,G,T,P}

    invprob = s.investmentproblem
    Rs = invprob.regionnames
    Gs = units.name
    Ts = string.(1:T)
    Ps = invprob.periodnames

    # Variables

    ops.committed = @variable(m, [1:R, 1:G, 1:T, 1:P], Int)
    varnames!(ops.committed, "committed_thermal_$(s.name)", Rs, Gs, Ts, Ps)

    ops.started   = @variable(m, [1:R, 1:G, 1:T, 1:P], Int)
    varnames!(ops.started, "started_thermal_$(s.name)", Rs, Gs, Ts, Ps)

    ops.shutdown  = @variable(m, [1:R, 1:G, 1:T, 1:P], Int)
    varnames!(ops.shutdown, "shutdown_thermal_$(s.name)", Rs, Gs, Ts, Ps)

    ops.energydispatch = @variable(m, [1:R, 1:G, 1:T, 1:P])
    varnames!(ops.energydispatch, "energydispatch_thermal_$(s.name)", Rs, Gs, Ts, Ps)

    ops.raisereserve   = @variable(m, [1:R, 1:G, 1:T, 1:P])
    varnames!(ops.raisereserve, "raisereserve_thermal_$(s.name)", Rs, Gs, Ts, Ps)

    ops.lowerreserve   = @variable(m, [1:R, 1:G, 1:T, 1:P])
    varnames!(ops.lowerreserve, "lowerreserve_thermal_$(s.name)", Rs, Gs, Ts, Ps)

    # Expressions

    ops.fixedcosts =
        @expression(m, [r in 1:R, g in 1:G],
                    ops.fixedcost[r,g] * invs.dispatching[r,g])

    ops.variablecosts =
        @expression(m, [r in 1:R, g in 1:G, p in 1:P],
                    sum(ops.variablecost[r,g] * ops.energydispatch[r,g,t,p] +
                        ops.startupcost[r,g] * ops.started[r,g,t,p] +
                        ops.shutdowncost[r,g] * ops.shutdown[r,g,t,p]
                        for t in 1:T))

    ops.operatingcosts =
        @expression(m, [r in 1:R, g in 1:G],
                    ops.fixedcosts[r,g] +
                    sum(ops.variablecosts[r,g,p] * periodweights[p] for p in 1:P))

    ops.ucap =
        @expression(m, [r in 1:R], G > 0 ? sum(
            invs.dispatching[r,g] * units.maxgen[g] * units.capacitycredit[g]
            for g in 1:G) : 0)

    ops.totalenergy =
        @expression(m, [r in 1:R, t in 1:T, p in 1:P],
                    G > 0 ? sum(ops.energydispatch[r,g,t,p] for g in 1:G) : 0)

    ops.totalraisereserve =
        @expression(m, [r in 1:R, t in 1:T, p in 1:P],
                    G > 0 ? sum(ops.raisereserve[r,g,t,p] for g in 1:G) : 0)

    ops.totallowerreserve =
        @expression(m, [r in 1:R, t in 1:T, p in 1:P],
                    G > 0 ? sum(ops.lowerreserve[r,g,t,p] for g in 1:G) : 0)

    # Constraints

    ops.minunitcommitment =
        @constraint(m, [r in 1:R, g in 1:G, t in 1:T, p in 1:P],
                    ops.committed[r,g,t,p] >= 0)

    ops.maxunitcommitment =
        @constraint(m, [r in 1:R, g in 1:G, t in 1:T, p in 1:P],
                    ops.committed[r,g,t,p] <= invs.dispatching[r,g])

    ops.minunitstartups =
        @constraint(m, [r in 1:R, g in 1:G, t in 1:T, p in 1:P],
                    ops.started[r,g,t,p] >= 0)

    ops.minunitshutdowns =
        @constraint(m, [r in 1:R, g in 1:G, t in 1:T, p in 1:P],
                    ops.shutdown[r,g,t,p] >= 0)

    ops.unitcommitmentcontinuity =
        @constraint(m, [r in 1:R, g in 1:G, t in 1:T, p in 1:P],
                    ops.committed[r,g,t,p] ==
                    ops.committed[r,g,prev(1,t,T),p]
                    + ops.started[r,g,t,p] - ops.shutdown[r,g,t,p])

    ops.minunituptime =
        @constraint(m, [r in 1:R, g in 1:G, t in 1:T, p in 1:P],
                    sum(ops.started[r,g,prev(i,t,T),p]
                        for i in 0:(units.minuptime[g]-1)) <=
                    ops.committed[r,g,t,p])

    ops.minunitdowntime =
        @constraint(m, [r in 1:R, g in 1:G, t in 1:T, p in 1:P],
                    sum(ops.shutdown[r,g,prev(i,t,T),p]
                        for i in 0:(units.mindowntime[g]-1)) <=
                    invs.dispatching[r,g] - ops.committed[r,g,t,p])

    ops.mingeneration =
        @constraint(m, [r in 1:R, g in 1:G, t in 1:T, p in 1:P],
                    ops.energydispatch[r,g,t,p] - ops.lowerreserve[r,g,t,p] >=
                    ops.committed[r,g,t,p] * units.mingen[g])

    ops.maxgeneration =
        @constraint(m, [r in 1:R, g in 1:G, t in 1:T, p in 1:P],
                    ops.energydispatch[r,g,t,p] + ops.raisereserve[r,g,t,p] <=
                    ops.committed[r,g,t,p] * units.maxgen[g])

    ops.minlowerreserve =
        @constraint(m, [r in 1:R, g in 1:G, t in 1:T, p in 1:P],
                    ops.lowerreserve[r,g,t,p] >= 0)

    ops.maxlowerreserve =
        @constraint(m, [r in 1:R, g in 1:G, t in 1:T, p in 1:P],
                    ops.lowerreserve[r,g,t,p] <=
                    ops.committed[r,g,t,p] * units.maxrampdown[g])

    ops.minraisereserve =
        @constraint(m, [r in 1:R, g in 1:G, t in 1:T, p in 1:P],
                    ops.raisereserve[r,g,t,p] >= 0)

    ops.maxraisereserve =
        @constraint(m, [r in 1:R, g in 1:G, t in 1:T, p in 1:P],
                    ops.raisereserve[r,g,t,p] <=
                    ops.committed[r,g,t,p] * units.maxrampup[g])

    ops.maxrampdown =
        @constraint(m, [r in 1:R, g in 1:G, t in 1:T, p in 1:P],
                    ops.energydispatch[r,g,t,p] - ops.lowerreserve[r,g,t,p] >=
                    ops.energydispatch[r,g,prev(1,t,T),p]
                    - (ops.committed[r,g,t,p] - ops.started[r,g,t,p]) *
                      units.maxrampdown[g]
                    + ops.started[r,g,t,p] * units.mingen[g]
                    - ops.shutdown[r,g,t,p] *
                      max(units.mingen[g], units.maxrampdown[g]))

    ops.maxrampup =
        @constraint(m, [r in 1:R, g in 1:G, t in 1:T, p in 1:P],
                    ops.energydispatch[r,g,t,p] + ops.lowerreserve[r,g,t,p] <=
                    ops.energydispatch[r,g,prev(1,t,T),p]
                    + (ops.committed[r,g,t,p] - ops.started[r,g,t,p]) *
                      units.maxrampup[g]
                    - ops.shutdown[r,g,t,p] * units.mingen[g]
                    + ops.started[r,g,t,p] *
                      max(units.mingen[g], units.maxrampup[g]))

end

welfare(x::ThermalGeneratorOperations) = -sum(x.operatingcosts)


text:
@testset "@fastcall" begin
    @test macroexpand(:(@fastcall f())) == :(f())
    @test macroexpand(:(@fastcall f(x))) == :(f(x))
    @test macroexpand(:(@fastcall f(x,y))) == :(f(x,y))

    # these are confusing
    @test_broken macroexpand(:(@fastcall f(a=1))) == :(GlobalRef(FastKeywords, Symbol("fastkw#f"))(GlobalRef(FastKeywords,:KW){:a}(1)))
    @test_broken macroexpand(:(@fastcall f(a=1,b=2))) == :(GlobalRef(FastKeywords, Symbol("fastkw#f"))(FastKeywords.KW{:a}(1), KW{:b}(2)))

end


text: largez01() = largez(1)


text:
using ITensors, Test

function is_unitary(U::ITensor; kwargs...)
  s = noprime(filterinds(U; plev=1))
  return isapprox(transpose(dag(U))(U), op("I", s...))
end

@testset "SiteType" begin
  N = 10

  @testset "Star in operator strings" begin
    sites = siteinds("S=1/2", N)
    #@test_throws ArgumentError op(sites, "Sp", 1)
    @test sites[1] isa Index
    Sz = op(sites, "Sz", 2)
    SzSz = op(sites, "Sz * Sz", 2)
    @test SzSz ≈ product(Sz, Sz)
    Sy = op(sites, "Sy", 2)
    SySy = op(sites, "Sy * Sy", 2)
    @test SySy ≈ product(Sy, Sy)

    Sz1 = op("Sz", sites, 1)
    @test op("Sz", [sites[1]]) ≈ Sz1
    @test op([sites[1]], "Sz") ≈ Sz1
    @test op([1 0; 0 -1] / 2, [sites[1]]) ≈ Sz1
    @test op([sites[1]], [1 0; 0 -1] / 2) ≈ Sz1

    @test op([sites[1]], "Ry"; θ=π / 2) ≈
      itensor([1 -1; 1 1] / √2, sites[1]', dag(sites[1]))

    sites = siteinds("S=1", N)
    #@test_throws ArgumentError op(sites, "Sp", 1)
    Sz = op(sites, "Sz", 2)
    SzSz = op(sites, "Sz * Sz", 2)
    @test SzSz ≈ product(Sz, Sz)
    Sy = op(sites, "Sy", 2)
    SySy = op(sites, "Sy * Sy", 2)
    @test SySy ≈ product(Sy, Sy)
    SzSySz = op(sites, "Sz * Sy * Sz", 2)
    @test SzSySz ≈ product(Sz, product(Sy, Sz))
  end

  @testset "+/- in operator strings" begin
    q = siteind("Qudit"; dim=5)
    Amat = array(op("a", q))
    Adagmat = array(op("a†", q))

    x = Amat - Adagmat
    @test x ≈ array(op("a - a†", q))
    x = Amat * Adagmat - Adagmat
    @test x ≈ array(op("a * a† - a†", q))
    @test x ≈ array(op("a * a† - a†", q))
    x = Adagmat * Adagmat * Amat * Amat
    @test x ≈ array(op("a† * a† * a * a", q))

    q = siteind("S=1/2")
    Sp = array(op("S+", q))
    Sm = array(op("S-", q))
    Sx = array(op("Sx", q))
    Sy = array(op("Sy", q))
    Sz = array(op("Sz", q))
    x = Sp + Sm
    @test x ≈ array(op("S+ + S-", q))
    x = Sp - Sm
    @test x ≈ array(op("S+ - S-", q))
    x = Sp - Sm - Sp
    @test x ≈ array(op("S+ - S- - S+", q))
    x = Sp * Sm + Sm * Sp
    @test x ≈ array(op("S+ * S- + S- * S+", q))
    # Deprecated syntax
    @test x ≈ array(op("S+ * S- + S-*S+", q))
    x = Sp * Sm - Sm * Sp
    @test x ≈ array(op("S+ * S- - S- * S+", q))
    @test x ≈ array(op("S+ * S- - S- * S+", q))
    x = Sp * Sm + Sm * Sp + Sz * Sx * Sy
    @test x ≈ array(op("S+ * S- + S- * S+ + Sz * Sx * Sy", q))
    x = Sp * Sm - Sm * Sp + Sz * Sx * Sy
    @test x ≈ array(op("S+ * S- - S- * S+ + Sz * Sx * Sy", q))
    x = Sp * Sm - Sm * Sp - Sz * Sx * Sy
    @test x ≈ array(op("S+ * S- - S- * S+ - Sz * Sx * Sy", q))

    #q = siteind("Qubit")
    #R = array(op("Rx", q; θ = 0.1))
    #H = array(op("H", q))
    #Y = array(op("Y", q))
    #x = H * R + Y + R
    #@test x ≈ array(op("H * Rx + Y + Rx", q; θ = 0.1))

  end

  @testset "Custom SiteType using op" begin
    # Use "_Custom_" tag even though this example
    # is for S=3/2, because we might define the
    # "S=3/2" TagType inside ITensors.jl later
    function ITensors.op(::OpName"Sz", ::SiteType"_Custom_", s::Index)
      Op = emptyITensor(s', dag(s))
      Op[s' => 1, s => 1] = +3 / 2
      Op[s' => 2, s => 2] = +1 / 2
      Op[s' => 3, s => 3] = -1 / 2
      Op[s' => 4, s => 4] = -3 / 2
      return Op
    end

    function ITensors.op(::OpName"α", ::SiteType"_Custom_", s1::Index, s2::Index)
      Op = emptyITensor(s1', s2', dag(s1), dag(s2))
      Op[s1' => 1, s2' => 2, s1 => 1, s2 => 2] = +3 / 2
      Op[s1' => 2, s2' => 1, s1 => 2, s2 => 2] = +1 / 2
      Op[s1' => 3, s2' => 3, s1 => 3, s2 => 4] = -1 / 2
      Op[s1' => 4, s2' => 1, s1 => 4, s2 => 2] = -3 / 2
      return Op
    end

    function ITensors.op(
      ::OpName"β", ::SiteType"_Custom1", ::SiteType"_Custom2", s1::Index, s2::Index
    )
      Op = emptyITensor(s1', s2', dag(s1), dag(s2))
      Op[s1' => 1, s2' => 2, s1 => 1, s2 => 2] = +5 / 2
      Op[s1' => 2, s2' => 1, s1 => 2, s2 => 2] = +3 / 2
      Op[s1' => 3, s2' => 3, s1 => 3, s2 => 4] = -3 / 2
      Op[s1' => 4, s2' => 1, s1 => 4, s2 => 2] = -5 / 2
      return Op
    end

    s = Index(4, "_Custom_, __x")
    Sz = op("Sz", s)
    @test Sz[s' => 1, s => 1] ≈ +3 / 2
    @test Sz[s' => 2, s => 2] ≈ +1 / 2
    @test Sz[s' => 3, s => 3] ≈ -1 / 2
    @test Sz[s' => 4, s => 4] ≈ -3 / 2

    t = Index(4, "_Custom_, __x")
    α = op("α", s, t)
    @test α[s' => 1, t' => 2, s => 1, t => 2] ≈ +3 / 2
    @test α[s' => 2, t' => 1, s => 2, t => 2] ≈ +1 / 2
    @test α[s' => 3, t' => 3, s => 3, t => 4] ≈ -1 / 2
    @test α[s' => 4, t' => 1, s => 4, t => 2] ≈ -3 / 2

    s1 = Index(4, "_Custom1, __x")
    @test_throws ArgumentError op("α", s, s1)

    s2 = Index(4, "_Custom2, __x")
    β = op("β", s1, s2)
    @test β[s1' => 1, s2' => 2, s1 => 1, s2 => 2] ≈ +5 / 2
    @test β[s1' => 2, s2' => 1, s1 => 2, s2 => 2] ≈ +3 / 2
    @test β[s1' => 3, s2' => 3, s1 => 3, s2 => 4] ≈ -3 / 2
    @test β[s1' => 4, s2' => 1, s1 => 4, s2 => 2] ≈ -5 / 2
    @test_throws ArgumentError op("β", s2, s1)
  end

  @testset "Custom OpName with long name" begin
    function ITensors.op(::OpName"my_favorite_operator", ::SiteType"S=1/2", s::Index)
      Op = emptyITensor(s', dag(s))
      Op[s' => 1, s => 1] = 0.11
      Op[s' => 1, s => 2] = 0.12
      Op[s' => 2, s => 1] = 0.21
      Op[s' => 2, s => 2] = 0.22
      return Op
    end

    s = Index(2, "S=1/2, Site")
    Sz = op("my_favorite_operator", s)
    @test Sz[s' => 1, s => 1] ≈ 0.11
    @test Sz[s' => 1, s => 2] ≈ 0.12
    @test Sz[s' => 2, s => 1] ≈ 0.21
    @test Sz[s' => 2, s => 2] ≈ 0.22
  end

  @testset "op with more than two indices" begin
    ITensors.space(::SiteType"qubit") = 2

    function ITensors.op(::OpName"rand", ::SiteType"qubit", s::Index...)
      return randomITensor(prime.(s)..., dag.(s)...)
    end

    s = siteinds("qubit", 4)
    o = op("rand", s...)
    @test norm(o) > 0
    @test order(o) == 8
    @test hassameinds(o, (prime.(s)..., s...))
  end

  @testset "Custom SiteType using op!" begin
    # Use "_Custom_" tag even though this example
    # is for S=3/2, because we might define the
    # "S=3/2" TagType inside ITensors.jl later
    function ITensors.op!(Op::ITensor, ::OpName"Sz", ::SiteType"_Custom_", s::Index)
      Op[s' => 1, s => 1] = +3 / 2
      Op[s' => 2, s => 2] = +1 / 2
      Op[s' => 3, s => 3] = -1 / 2
      return Op[s' => 4, s => 4] = -3 / 2
    end

    function ITensors.op!(
      Op::ITensor, ::OpName"α", ::SiteType"_Custom_", s1::Index, s2::Index
    )
      Op[s1' => 1, s2' => 2, s1 => 1, s2 => 2] = +3 / 2
      Op[s1' => 2, s2' => 1, s1 => 2, s2 => 2] = +1 / 2
      Op[s1' => 3, s2' => 3, s1 => 3, s2 => 4] = -1 / 2
      return Op[s1' => 4, s2' => 1, s1 => 4, s2 => 2] = -3 / 2
    end

    function ITensors.op!(
      Op::ITensor,
      ::OpName"β",
      ::SiteType"_Custom1",
      ::SiteType"_Custom2",
      s1::Index,
      s2::Index,
    )
      Op[s1' => 1, s2' => 2, s1 => 1, s2 => 2] = +5 / 2
      Op[s1' => 2, s2' => 1, s1 => 2, s2 => 2] = +3 / 2
      Op[s1' => 3, s2' => 3, s1 => 3, s2 => 4] = -3 / 2
      return Op[s1' => 4, s2' => 1, s1 => 4, s2 => 2] = -5 / 2
    end

    s = Index(4, "_Custom_, __x")
    Sz = op("Sz", s)
    @test Sz[s' => 1, s => 1] ≈ +3 / 2
    @test Sz[s' => 2, s => 2] ≈ +1 / 2
    @test Sz[s' => 3, s => 3] ≈ -1 / 2
    @test Sz[s' => 4, s => 4] ≈ -3 / 2

    t = Index(4, "_Custom_, __x")
    α = op("α", s, t)
    @test α[s' => 1, t' => 2, s => 1, t => 2] ≈ +3 / 2
    @test α[s' => 2, t' => 1, s => 2, t => 2] ≈ +1 / 2
    @test α[s' => 3, t' => 3, s => 3, t => 4] ≈ -1 / 2
    @test α[s' => 4, t' => 1, s => 4, t => 2] ≈ -3 / 2

    s1 = Index(4, "_Custom1, __x")
    @test_throws ArgumentError op("α", t, s1)

    s2 = Index(4, "_Custom2, __x")
    β = op("β", s1, s2)
    @test β[s1' => 1, s2' => 2, s1 => 1, s2 => 2] ≈ +5 / 2
    @test β[s1' => 2, s2' => 1, s1 => 2, s2 => 2] ≈ +3 / 2
    @test β[s1' => 3, s2' => 3, s1 => 3, s2 => 4] ≈ -3 / 2
    @test β[s1' => 4, s2' => 1, s1 => 4, s2 => 2] ≈ -5 / 2
    @test_throws ArgumentError op("β", s2, s1)
  end

  @testset "Custom SiteType using older op interface" begin
    # Use "_Custom_" tag even though this example
    # is for S=3/2, because we might define the
    # "S=3/2" TagType inside ITensors.jl later
    function ITensors.op(::SiteType"_Custom_", s::Index, opname::AbstractString)
      Op = emptyITensor(s', dag(s))
      if opname == "S+"
        Op[s' => 1, s => 2] = sqrt(3)
        Op[s' => 2, s => 3] = 2
        Op[s' => 3, s => 4] = sqrt(3)
      else
        error("Name $opname not recognized for tag \"Custom\"")
      end
      return Op
    end

    s = Index(4, "_Custom_")
    Sp = op("S+", s)
    @test Sp[s' => 1, s => 2] ≈ sqrt(3)
    @test Sp[s' => 2, s => 3] ≈ 2
    @test Sp[s' => 3, s => 4] ≈ sqrt(3)
  end

  @testset "siteind defined by space overload" begin
    ITensors.space(::SiteType"Test1") = 4
    s = siteind("Test1", 3)
    @test dim(s) == 4
    @test hastags(s, "Site,Test1,n=3")

    s = siteind("Test1")
    @test dim(s) == 4
    @test hastags(s, "Site,Test1")
  end

  @testset "siteind defined by siteind overload" begin
    ITensors.siteind(::SiteType"Test2") = Index(4, "Test2")
    s = siteind("Test2", 3)
    @test dim(s) == 4
    @test hastags(s, "Test2,n=3")
  end

  @testset "siteind defined by space overload with QN" begin
    function ITensors.space(::SiteType"Test3")
      return [QN("T", 0) => 2, QN("T", 1) => 1, QN("T", 2) => 1]
    end
    s = siteind("Test3", 3)
    @test dim(s) == 4
    @test hasqns(s)
    @test hastags(s, "Site,Test3,n=3")
  end

  @testset "siteinds defined by space overload" begin
    function ITensors.space(::SiteType"Test4"; conserve_qns=false)
      if conserve_qns
        return [QN("T", 0) => 2, QN("T", 1) => 1, QN("T", 2) => 1]
      end
      return 4
    end

    # Without QNs
    s = siteinds("Test4", 6)
    @test length(s) == 6
    @test dim(s[1]) == 4
    for n in 1:length(s)
      @test hastags(s[n], "Site,Test4,n=$n")
      @test !hasqns(s[n])
    end

    # With QNs
    s = siteinds("Test4", 6; conserve_qns=true)
    @test length(s) == 6
    @test dim(s[1]) == 4
    for n in 1:length(s)
      @test hastags(s[n], "Site,Test4,n=$n")
      @test hasqns(s[n])
    end
  end

  @testset "siteinds defined by siteinds overload" begin
    function ITensors.siteinds(::SiteType"Test5", N; kwargs...)
      return [Index(4, "Test5,n=$n") for n in 1:N]
    end
    s = siteinds("Test5", 8)
    @test length(s) == 8
    @test dim(s[1]) == 4
    for n in 1:length(s)
      @test hastags(s[n], "Test5,n=$n")
    end
  end

  @testset "Version of siteinds taking function argument" begin
    N = 10
    s = siteinds(n -> (n == 1 || n == N) ? "S=1/2" : "S=1", N)
    for n in (1, N)
      @test dim(s[n]) == 2
      @test hastags(s[n], "Site,S=1/2,n=$n")
    end
    for n in 2:(N - 1)
      @test dim(s[n]) == 3
      @test hastags(s[n], "Site,S=1,n=$n")
    end
  end

  @testset "siteinds addtags keyword argument" begin
    N = 4
    s = siteinds("S=1/2", N; addtags="T")
    for n in 1:N
      @test hastags(s[n], "Site,S=1/2,n=$n,T")
    end
  end

  @testset "Error for undefined tag in siteinds,space system" begin
    @test_throws ErrorException siteinds("Missing", 10)
    @test_throws ErrorException siteind("Missing", 3)
    @test isnothing(siteind("Missing"))
  end

  @testset "Various ops input types" begin
    s = siteinds("S=1/2", 4)

    # Vector{Tuple{String,Int}} input
    oa = ops(s, [("Sz", n) for n in 1:length(s)])
    @test length(oa) == length(s)
    @test norm(oa[2] - op("Sz", s, 2)) < 1E-8

    # Vector{Tuple} input
    oa = ops(s, Tuple[("Sz", n) for n in 1:length(s)])
    @test length(oa) == length(s)
    @test norm(oa[2] - op("Sz", s, 2)) < 1E-8
  end

  @testset "Index Values From Strings" begin
    @testset "Val function" begin
      s = siteind("Electron")
      @test val(s, "0") == 1
      @test val(s, "Up") == 2
      @test val(s, "Dn") == 3
      @test val(s, "UpDn") == 4
    end

    @testset "Strings in ITensor get and set" begin
      s = siteind("S=1"; conserve_qns=true)
      T = ITensor(s', dag(s))
      T[s' => "Up", s => "Up"] = +1.0
      T[s' => "Z0", s => "Z0"] = +2.0
      T[s' => "Dn", s => "Dn"] = -1.0
      @test T[1, 1] ≈ +1.0
      @test T[2, 2] ≈ +2.0
      @test T[3, 3] ≈ -1.0

      o = onehot(s => "Z0")
      @test vector(o) ≈ [0, 1, 0]
    end
  end

  @testset "state with variable dimension" begin
    ITensors.space(::SiteType"MyQudit"; dim=2) = dim

    function ITensors.state(::StateName{N}, ::SiteType"MyQudit", s::Index) where {N}
      n = parse(Int, String(N))
      st = zeros(dim(s))
      st[n + 1] = 1.0
      return itensor(st, s)
    end

    s = siteind("MyQudit"; dim=3)
    v0 = state(s, "0")
    v1 = state(s, "1")
    v2 = state(s, "2")
    @test v0 == state("0", s)
    @test v1 == state("1", s)
    @test v2 == state("2", s)
    @test dim(v0) == 3
    @test dim(v1) == 3
    @test dim(v2) == 3
    @test v0[s => 1] == 1
    @test v0[s => 2] == 0
    @test v0[s => 3] == 0
    @test v1[s => 1] == 0
    @test v1[s => 2] == 1
    @test v1[s => 3] == 0
    @test v2[s => 1] == 0
    @test v2[s => 2] == 0
    @test v2[s => 3] == 1
    @test_throws BoundsError state(s, "3")
  end

  @testset "state with variable dimension (deprecated)" begin
    ITensors.space(::SiteType"MyQudit2"; dim=2) = dim

    # XXX: This syntax is deprecated, only testing for
    # backwards compatibility. Should return the
    # ITensor `itensor(st, s)`.
    function ITensors.state(::StateName{N}, ::SiteType"MyQudit2", s::Index) where {N}
      n = parse(Int, String(N))
      st = zeros(dim(s))
      st[n + 1] = 1.0
      return st
    end

    s = siteind("MyQudit2"; dim=3)
    v0 = state(s, "0")
    v1 = state(s, "1")
    v2 = state(s, "2")
    @test v0 == state("0", s)
    @test v1 == state("1", s)
    @test v2 == state("2", s)
    @test dim(v0) == 3
    @test dim(v1) == 3
    @test dim(v2) == 3
    @test v0[s => 1] == 1
    @test v0[s => 2] == 0
    @test v0[s => 3] == 0
    @test v1[s => 1] == 0
    @test v1[s => 2] == 1
    @test v1[s => 3] == 0
    @test v2[s => 1] == 0
    @test v2[s => 2] == 0
    @test v2[s => 3] == 1
    @test_throws BoundsError state(s, "3")
  end

  @testset "Regression test for state overload" begin
    ITensors.space(::SiteType"Xev") = 8
    function ITensors.state(::StateName"0", ::SiteType"Xev")
      return [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    end
    s = siteind("Xev")
    @test state(s, "0") ≈ ITensor([1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], s)
  end

  @testset "function applied to a gate" begin
    s = siteinds("Qubit", 2)

    θ = 0.1
    rx = array(op("Rx", s[1]; θ=0.1))
    exp_rx = exp(rx)
    gtest = op(x -> exp(x), "Rx", s[1]; θ=0.1)
    @test exp_rx ≈ array(op(x -> exp(x), "Rx", s[1]; θ=0.1))
    @test exp_rx ≈ array(op(x -> exp(x), ("Rx", 1, (θ=0.1,)), s))

    cx = 0.1 * reshape(array(op("CX", s[1], s[2])), (4, 4))
    exp_cx = reshape(exp(cx), (2, 2, 2, 2))
    @test exp_cx ≈ array(op(x -> exp(0.1 * x), "CX", s[1], s[2]))
    @test exp_cx ≈ array(op(x -> exp(0.1 * x), ("CX", (1, 2)), s))
  end

  @testset "Haar-random unitary RandomUnitary" begin
    s = siteinds(2, 3)

    U = op("RandomUnitary", s, 1, 2)
    @test eltype(U) == ComplexF64
    @test order(U) == 4
    @test is_unitary(U; rtol=1e-15)

    U = op("RandomUnitary", s, 1, 2, 3)
    @test eltype(U) == ComplexF64
    @test order(U) == 6
    @test is_unitary(U; rtol=1e-15)

    U = op("RandomUnitary", s, 1, 2; eltype=Float64)
    @test eltype(U) == Float64
    @test order(U) == 4
    @test is_unitary(U; rtol=1e-15)

    U = op("RandomUnitary", s, 1, 2, 3; eltype=Float64)
    @test eltype(U) == Float64
    @test order(U) == 6
    @test is_unitary(U; rtol=1e-15)
  end
end

nothing


text:
#    Copyright (c) 2021 J.A. Duffek
#    Copyright (c) 2000 D.M. Spink
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, see <http://www.gnu.org/licenses/>.

"""
    findspan(n::I,p::I,u::F,U::Vector{F})::I where {I<:Integer,F<:AbstractFloat}

Find the span of a B-Spline knot vector at a parametric point

# Arguments:
- `n`: nuber of control points `-1`
- `p`: spline degree
- `n`: parametric point
- `U`: knot sequence

# Output:
- `s`: knot span of `u` in `U`

# Examples:
```julia
julia> s = findspan(n,p,u,U)
```

# Reference:
- Algorithm A2.1 from 'The NURBS BOOK' pg68 ([1997Piegl](@cite))
"""
function findspan(n::I,p::I,u::F,U::Vector{F}
                 )::I where {I<:Integer,F<:AbstractFloat}
# TODO: the original C function passes the knot vector as reference

# special case
if u == U[n+2]; return n; end

# do binary search
low = p;
high = n + 1;
mid = (low + high) ÷ 2;
while (u < U[mid+1] || u >= U[mid+2])
  if u < U[mid+1]
    high = mid;
  else
    low = mid;
  end # if
  mid = (low + high) ÷ 2;
end # while
return mid;
end # findspan


text:
using Ion

suites = Ion.find_test_suites(Returns(true), pkgdir(Ion, "test", "examples"))

suites_proc = Ion.divide_tests(suites, 4)
suites_proc = Ion.divide_tests(suites, 3)

Ion.test(pkgdir(Ion, "test", "examples"), String["folder/*"])
Ion.test(pkgdir(Ion, "test", "examples"), String[])


text:
#ДАНО: Робот - в юго-западном углу поля, на котором расставлено некоторое количество маркеров
#РЕЗУЛЬТАТ: Функция вернула значение средней температуры всех замаркированных клеток

function check_temp(r) #функция возвращает температуру клетки + считает количество клеток с маркером
    global K
    if ismarker(r) == true
        K += 1
        return temperature(r)
    else
        return 0
    end
end

function sum_temp(r, side) #функция возвращает значение суммы температур всех клеток в направлении
    sum_temp = 0
    while isborder(r,side) == false
        sum_temp += check_temp(r)
        move!(r,side)
    end
    sum_temp += check_temp(r)
    return sum_temp
end

function reverse(side)
    for i=0:3
        if side == HorizonSide(i)
            return HorizonSide(mod(i+2, 4))
        end
    end
end

K = 0
function temp_corner(r)
    global K = 0
    temp = 0
    side = Ost
    while isborder(r,Nord) == false
        temp += sum_temp(r,side)
        move!(r,Nord)
        side = reverse(side)
    end
    temp += sum_temp(r,side)
    average_temp = temp/K
    return average_temp
end

max_stars_repo_name: gardend/mirea-progs

max_issues_repo_name: gardend/mirea-progs

max_forks_repo_name: gardend/mirea-progs