text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

#
# package typedb.client.test.behaviour.graql;
#
# import typedb.client.api.answer.ConceptMap;
# import typedb.client.api.answer.ConceptMapGroup;
# import typedb.client.api.answer.Numeric;
# import typedb.client.api.answer.NumericGroup;
# import typedb.client.api.concept.Concept;
# import typedb.client.api.concept.thing.Attribute;
# import typedb.client.api.concept.type.AttributeType;
# import typedb.client.api.query.QueryFuture;
# import typedb.client.common.Label;
# import graql.lang.Graql;
# import graql.lang.common.exception.GraqlException;
# import graql.lang.query.GraqlDefine;
# import graql.lang.query.GraqlDelete;
# import graql.lang.query.GraqlInsert;
# import graql.lang.query.GraqlMatch;
# import graql.lang.query.GraqlUndefine;
# import graql.lang.query.GraqlUpdate;
# import io.cucumber.java.en.Given;
# import io.cucumber.java.en.Then;
# import io.cucumber.java.en.When;
#
# import java.time.LocalDate;
# import java.time.LocalDateTime;
# import java.time.format.DateTimeParseException;
# import java.util.ArrayList;
# import java.util.HashMap;
# import java.util.List;
# import java.util.Map;
# import java.util.Set;
# import java.util.regex.Matcher;
# import java.util.regex.Pattern;
# import java.util.stream.Collectors;
# import java.util.stream.Stream;
#
# import static typedb.client.test.behaviour.connection.ConnectionStepsBase.tx;
# import static typedb.client.test.behaviour.util.Util.assertThrows;
# import static typedb.client.test.behaviour.util.Util.assertThrowsWithMessage;
# import static org.junit.Assert.assertEquals;
# import static org.junit.Assert.assertNotNull;
# import static org.junit.Assert.assertTrue;
#
# public class GraqlSteps {
#
#     private static List<ConceptMap> answers;
#     private static Numeric numericAnswer;
#     private static List<ConceptMapGroup> answerGroups;
#     private static List<NumericGroup> numericAnswerGroups;
#     private Map<String, Map<String, String>> rules;
#
#     @Given("the integrity is validated")
#     public void integrity_is_validated() {
#
#         // TODO
#
#     }
#
#     @Given("graql define")
#     public QueryFuture<Void> graql_define(String defineQueryStatements) {
#         GraqlDefine graqlQuery = Graql.parseQuery(String.join("\n", defineQueryStatements));
#         return tx().query().define(graqlQuery);
#     }
#
#     @Given("graql define; throws exception")
#     public void graql_define_throws(String defineQueryStatements) {
#         assertThrows(() -> graql_define(defineQueryStatements).get());
#     }
#
#     @Given("graql define; throws exception containing {string}")
#     public void graql_define_throws_exception(String exception, String defineQueryStatements) {
#         assertThrowsWithMessage(() -> graql_define(defineQueryStatements).get(), exception);
#     }
#
#     @Given("graql undefine")
#     public QueryFuture<Void> graql_undefine(String undefineQueryStatements) {
#         GraqlUndefine graqlQuery = Graql.parseQuery(String.join("\n", undefineQueryStatements));
#         return tx().query().undefine(graqlQuery);
#     }
#
#     @Given("graql undefine; throws exception")
#     public void graql_undefine_throws(String undefineQueryStatements) {
#         assertThrows(() -> graql_undefine(undefineQueryStatements).get());
#     }
#
#     @Given("graql undefine; throws exception containing {string}")
#     public void graql_undefine_throws_exception(String exception, String undefineQueryStatements) {
#         assertThrowsWithMessage(() -> graql_undefine(undefineQueryStatements).get(), exception);
#     }
#
#     @Given("graql insert")
#     public Stream<ConceptMap> graql_insert(String insertQueryStatements) {
#         GraqlInsert graqlQuery = Graql.parseQuery(String.join("\n", insertQueryStatements));
#         return tx().query().insert(graqlQuery);
#     }
#
#     @Given("graql insert; throws exception")
#     @SuppressWarnings("ReturnValueIgnored")
#     public void graql_insert_throws(String insertQueryStatements) {
#         //noinspection ResultOfMethodCallIgnored
#         assertThrows(() -> graql_insert(insertQueryStatements).collect(Collectors.toList()));
#     }
#
#     @Given("graql insert; throws exception containing {string}")
#     @SuppressWarnings("ReturnValueIgnored")
#     public void graql_insert_throws_exception(String exception, String insertQueryStatements) {
#         //noinspection ResultOfMethodCallIgnored
#         assertThrowsWithMessage(() -> graql_insert(insertQueryStatements).collect(Collectors.toList()), exception);
#     }
#
#     @Given("graql delete")
#     public QueryFuture<Void> graql_delete(String deleteQueryStatements) {
#         GraqlDelete graqlQuery = Graql.parseQuery(String.join("\n", deleteQueryStatements));
#         return tx().query().delete(graqlQuery);
#     }
#
#     @Given("graql delete; throws exception")
#     public void graql_delete_throws(String deleteQueryStatements) {
#         assertThrows(() -> graql_delete(deleteQueryStatements).get());
#     }
#
#     @Given("graql delete; throws exception containing {string}")
#     public void graql_delete_throws_exception(String exception, String deleteQueryStatements) {
#         assertThrowsWithMessage(() -> graql_delete(deleteQueryStatements).get(), exception);
#     }
#
#     @Given("graql update")
#     public Stream<ConceptMap> graql_update(String updateQueryStatements) {
#         GraqlUpdate graqlQuery = Graql.parseQuery(String.join("\n", updateQueryStatements));
#         return tx().query().update(graqlQuery);
#     }
#
#     @Given("graql update; throws exception")
#     public void graql_update_throws(String updateQueryStatements) {
#         assertThrows(() -> graql_update(updateQueryStatements).iterator().next());
#     }
#
#     @Given("graql update; throws exception containing {string}")
#     public void graql_update_throws_exception(String exception, String updateQueryStatements) {
#         assertThrowsWithMessage(() -> graql_update(updateQueryStatements).iterator().next(), exception);
#     }
#
#     private void clearAnswers() {
#         answers = null;
#         numericAnswer = null;
#         answerGroups = null;
#         numericAnswerGroups = null;
#     }
#
#     @When("get answers of graql insert")
#     public void get_answers_of_graql_insert(String graqlQueryStatements) {
#         GraqlInsert graqlQuery = Graql.parseQuery(String.join("\n", graqlQueryStatements));
#         clearAnswers();
#         answers = tx().query().insert(graqlQuery).collect(Collectors.toList());
#     }
#
#     @When("get answers of graql match")
#     public void graql_match(String graqlQueryStatements) {
#         try {
#             GraqlMatch graqlQuery = Graql.parseQuery(String.join("\n", graqlQueryStatements)).asMatch();
#             clearAnswers();
#             answers = tx().query().match(graqlQuery).collect(Collectors.toList());
#         } catch (GraqlException e) {
#             // NOTE: We manually close transaction here, because we want to align with all non-java clients,
#             // where parsing happens at server-side which closes transaction if they fail
#             tx().close();
#             throw e;
#         }
#     }
#
#     @When("graql match; throws exception")
#     public void graql_match_throws_exception(String graqlQueryStatements) {
#         assertThrows(() -> graql_match(graqlQueryStatements));
#     }
#
#     @When("get answer of graql match aggregate")
#     public void graql_match_aggregate(String graqlQueryStatements) {
#         GraqlMatch.Aggregate graqlQuery = Graql.parseQuery(String.join("\n", graqlQueryStatements)).asMatchAggregate();
#         clearAnswers();
#         numericAnswer = tx().query().match(graqlQuery).get();
#     }
#
#     @When("get answers of graql match group")
#     public void graql_match_group(String graqlQueryStatements) {
#         GraqlMatch.Group graqlQuery = Graql.parseQuery(String.join("\n", graqlQueryStatements)).asMatchGroup();
#         clearAnswers();
#         answerGroups = tx().query().match(graqlQuery).collect(Collectors.toList());
#     }
#
#     @When("get answers of graql match group aggregate")
#     public void graql_match_group_aggregate(String graqlQueryStatements) {
#         GraqlMatch.Group.Aggregate graqlQuery = Graql.parseQuery(String.join("\n", graqlQueryStatements)).asMatchGroupAggregate();
#         clearAnswers();
#         numericAnswerGroups = tx().query().match(graqlQuery).collect(Collectors.toList());
#     }
#
#     @Then("answer size is: {number}")
#     public void answer_quantity_assertion(int expectedAnswers) {
#         assertEquals(String.format("Expected [%d] answers, but got [%d]", expectedAnswers, answers.size()),
#                      expectedAnswers, answers.size());
#     }
#
#     @Then("uniquely identify answer concepts")
#     public void uniquely_identify_answer_concepts(List<Map<String, String>> answerConcepts) {
#         assertEquals(
#                 String.format("The number of identifier entries (rows) should match the number of answers, but found %d identifier entries and %d answers.",
#                               answerConcepts.size(), answers.size()),
#                 answerConcepts.size(), answers.size()
#         );
#
#         for (ConceptMap answer : answers) {
#             List<Map<String, String>> matchingIdentifiers = new ArrayList<>();
#
#             for (Map<String, String> answerIdentifier : answerConcepts) {
#
#                 if (matchAnswerConcept(answerIdentifier, answer)) {
#                     matchingIdentifiers.add(answerIdentifier);
#                 }
#             }
#             assertEquals(
#                     String.format("An identifier entry (row) should match 1-to-1 to an answer, but there were %d matching identifier entries for answer with variables %s.",
#                                   matchingIdentifiers.size(), answer.map().keySet().toString()),
#                     1, matchingIdentifiers.size()
#             );
#         }
#     }
#
#     @Then("order of answer concepts is")
#     public void order_of_answer_concepts_is(List<Map<String, String>> answersIdentifiers) {
#         assertEquals(
#                 String.format("The number of identifier entries (rows) should match the number of answers, but found %d identifier entries and %d answers.",
#                               answersIdentifiers.size(), answers.size()),
#                 answersIdentifiers.size(), answers.size()
#         );
#         for (int i = 0; i < answers.size(); i++) {
#             ConceptMap answer = answers.get(i);
#             Map<String, String> answerIdentifiers = answersIdentifiers.get(i);
#             assertTrue(
#                     String.format("The answer at index %d does not match the identifier entry (row) at index %d.", i, i),
#                     matchAnswerConcept(answerIdentifiers, answer)
#             );
#         }
#     }
#
#     @Then("aggregate value is: {double}")
#     public void aggregate_value_is(double expectedAnswer) {
#         assertNotNull("The last executed query was not an aggregate query", numericAnswer);
#         assertEquals(String.format("Expected answer to equal %f, but it was %f.", expectedAnswer, numericAnswer.asNumber().doubleValue()),
#                      expectedAnswer, numericAnswer.asNumber().doubleValue(), 0.001);
#     }
#
#     @Then("aggregate answer is not a number")
#     public void aggregate_answer_is_not_a_number() {
#         assertTrue(numericAnswer.isNaN());
#     }
#
#     @Then("answer groups are")
#     public void answer_groups_are(List<Map<String, String>> answerIdentifierTable) {
#         Set<AnswerIdentifierGroup> answerIdentifierGroups = answerIdentifierTable.stream()
#                 .collect(Collectors.groupingBy(x -> x.get(AnswerIdentifierGroup.GROUP_COLUMN_NAME)))
#                 .values()
#                 .stream()
#                 .map(AnswerIdentifierGroup::new)
#                 .collect(Collectors.toSet());
#
#         assertEquals(String.format("Expected [%d] answer groups, but found [%d].",
#                                    answerIdentifierGroups.size(), answerGroups.size()),
#                      answerIdentifierGroups.size(), answerGroups.size()
#         );
#
#         for (AnswerIdentifierGroup answerIdentifierGroup : answerIdentifierGroups) {
#             String[] identifier = answerIdentifierGroup.ownerIdentifier.split(":", 2);
#             UniquenessCheck checker;
#             switch (identifier[0]) {
#                 case "label":
#                     checker = new LabelUniquenessCheck(identifier[1]);
#                     break;
#                 case "key":
#                     checker = new KeyUniquenessCheck(identifier[1]);
#                     break;
#                 case "value":
#                     checker = new ValueUniquenessCheck(identifier[1]);
#                     break;
#                 default:
#                     throw new IllegalStateException("Unexpected value: " + identifier[0]);
#             }
#             ConceptMapGroup answerGroup = answerGroups.stream()
#                     .filter(ag -> checker.check(ag.owner()))
#                     .findAny()
#                     .orElse(null);
#             assertNotNull(String.format("The group identifier [%s] does not match any of the answer group owners.", answerIdentifierGroup.ownerIdentifier), answerGroup);
#
#             List<Map<String, String>> answersIdentifiers = answerIdentifierGroup.answersIdentifiers;
#             for (ConceptMap answer : answerGroup.conceptMaps()) {
#                 List<Map<String, String>> matchingIdentifiers = new ArrayList<>();
#
#                 for (Map<String, String> answerIdentifiers : answersIdentifiers) {
#
#                     if (matchAnswerConcept(answerIdentifiers, answer)) {
#                         matchingIdentifiers.add(answerIdentifiers);
#                     }
#                 }
#                 assertEquals(
#                         String.format("An identifier entry (row) should match 1-to-1 to an answer, but there were [%d] matching identifier entries for answer with variables %s.",
#                                       matchingIdentifiers.size(), answer.map().keySet().toString()),
#                         1, matchingIdentifiers.size()
#                 );
#             }
#         }
#     }
#
#     @Then("group aggregate values are")
#     public void group_aggregate_values_are(List<Map<String, String>> answerIdentifierTable) {
#         Map<String, Double> expectations = new HashMap<>();
#         for (Map<String, String> answerIdentifierRow : answerIdentifierTable) {
#             String groupOwnerIdentifier = answerIdentifierRow.get(AnswerIdentifierGroup.GROUP_COLUMN_NAME);
#             double expectedAnswer = Double.parseDouble(answerIdentifierRow.get("value"));
#             expectations.put(groupOwnerIdentifier, expectedAnswer);
#         }
#
#         assertEquals(String.format("Expected [%d] answer groups, but found [%d].", expectations.size(), numericAnswerGroups.size()),
#                      expectations.size(), numericAnswerGroups.size()
#         );
#
#         for (Map.Entry<String, Double> expectation : expectations.entrySet()) {
#             String[] identifier = expectation.getKey().split(":", 2);
#             UniquenessCheck checker;
#             switch (identifier[0]) {
#                 case "label":
#                     checker = new LabelUniquenessCheck(identifier[1]);
#                     break;
#                 case "key":
#                     checker = new KeyUniquenessCheck(identifier[1]);
#                     break;
#                 case "value":
#                     checker = new ValueUniquenessCheck(identifier[1]);
#                     break;
#                 default:
#                     throw new IllegalStateException("Unexpected value: " + identifier[0]);
#             }
#             double expectedAnswer = expectation.getValue();
#             NumericGroup answerGroup = numericAnswerGroups.stream()
#                     .filter(ag -> checker.check(ag.owner()))
#                     .findAny()
#                     .orElse(null);
#             assertNotNull(String.format("The group identifier [%s] does not match any of the answer group owners.", expectation.getKey()), answerGroup);
#
#             double actualAnswer = answerGroup.numeric().asNumber().doubleValue();
#             assertEquals(
#                     String.format("Expected answer [%f] for group [%s], but got [%f]",
#                                   expectedAnswer, expectation.getKey(), actualAnswer),
#                     expectedAnswer, actualAnswer, 0.001
#             );
#         }
#     }
#
#     @Then("number of groups is: {int}")
#     public void number_of_groups_is(int expectedGroupCount) {
#         assertEquals(expectedGroupCount, answerGroups.size());
#     }
#
#     public static class AnswerIdentifierGroup {
#         private final String ownerIdentifier;
#         private final List<Map<String, String>> answersIdentifiers;
#
#         private static final String GROUP_COLUMN_NAME = "owner";
#
#         public AnswerIdentifierGroup(List<Map<String, String>> answerIdentifierTable) {
#             ownerIdentifier = answerIdentifierTable.get(0).get(GROUP_COLUMN_NAME);
#             answersIdentifiers = new ArrayList<>();
#             for (Map<String, String> rawAnswerIdentifiers : answerIdentifierTable) {
#                 answersIdentifiers.add(rawAnswerIdentifiers.entrySet().stream()
#                                                .filter(e -> !e.getKey().equals(GROUP_COLUMN_NAME))
#                                                .collect(Collectors.toMap(Map.Entry::getKey, Map.Entry::getValue)));
#             }
#         }
#     }
#
#     // This method is used by explanation steps
# //    private boolean matchAnswer(Map<String, String> answerIdentifiers, ConceptMap answer) {
# //
# //        for (Map.Entry<String, String> entry : answerIdentifiers.entrySet()) {
# //            final Reference.Name var = Reference.named(entry.getKey());
# //            final String identifier = entry.getValue();
# //
# //            if (!identifierChecks.containsKey(identifier)) {
# //                throw new ScenarioDefinitionException(String.format("Identifier \"%s\" hasn't previously been declared.", identifier));
# //            }
# //
# //            if (!identifierChecks.get(identifier).check(answer.get(var))) {
# //                return false;
# //            }
# //        }
# //        return true;
# //    }
#
#     private boolean matchAnswerConcept(Map<String, String> answerIdentifiers, ConceptMap answer) {
#         for (Map.Entry<String, String> entry : answerIdentifiers.entrySet()) {
#             String var = entry.getKey();
#             String[] identifier = entry.getValue().split(":", 2);
#             switch (identifier[0]) {
#                 case "label":
#                     if (!new LabelUniquenessCheck(identifier[1]).check(answer.get(var))) {
#                         return false;
#                     }
#                     break;
#                 case "key":
#                     if (!new KeyUniquenessCheck(identifier[1]).check(answer.get(var))) {
#                         return false;
#                     }
#                     break;
#                 case "value":
#                     if (!new ValueUniquenessCheck(identifier[1]).check(answer.get(var))) {
#                         return false;
#                     }
#                     break;
#             }
#         }
#         return true;
#     }
#
#     @Then("rules are")
#     public void rules_are(Map<String, Map<String, String>> rules) {
#         this.rules = rules;
#     }
#
#     @Then("rules contain: {type_label}")
#     public void rules_contain(String ruleLabel) {
#         assert (tx().logic().getRules().anyMatch(rule -> rule.getLabel().equals(ruleLabel)));
#     }
#
#     @Then("rules do not contain: {type_label}")
#     public void rules_do_not_contain(String ruleLabel) {
#         assert (tx().logic().getRules().noneMatch(rule -> rule.getLabel().equals(ruleLabel)));
#     }
#
#     @Then("answers contain explanation tree")
#     public void answers_contain_explanation_tree(Map<Integer, Map<String, String>> explanationTree) {
#         // TODO
#         throw new UnsupportedOperationException();
# //        checkExplanationEntry(answers, explanationTree, 0);
#     }
#
#     /* private void checkExplanationEntry(List<ConceptMap> answers, Map<Integer, Map<String, String>> explanationTree, Integer entryId) {
#         Map<String, String> explanationEntry = explanationTree.get(entryId);
#         String[] vars = explanationEntry.get("vars").split(", ");
#         String[] identifiers = explanationEntry.get("identifiers").split(", ");
#         String[] children = explanationEntry.get("children").split(", ");
#
#         if (vars.length != identifiers.length) {
#             throw new ScenarioDefinitionException(String.format("vars and identifiers do not correspond for explanation entry %d. Found %d vars and %s identifiers", entryId, vars.length, identifiers.length));
#         }
#
#         Map<String, String> answerIdentifiers = IntStream.range(0, vars.length).boxed().collect(Collectors.toMap(i -> vars[i], i -> identifiers[i]));
#
#         Optional<ConceptMap> matchingAnswer = answers.stream().filter(answer -> matchAnswer(answerIdentifiers, answer)).findFirst();
#
#         assertTrue(String.format("No answer found for explanation entry %d that satisfies the vars and identifiers given", entryId), matchingAnswer.isPresent());
#         ConceptMap answer = matchingAnswer.get();
#
#         String queryWithIds = applyQueryTemplate(explanationEntry.get("pattern"), answer);
#         Conjunction<?> queryWithIdsConj = Graql.and(Graql.parsePatternList(queryWithIds));
#         assertEquals(
#                 String.format("Explanation entry %d has an incorrect pattern.\nExpected: %s\nActual: %s", entryId, queryWithIdsConj, answer.queryPattern()),
#                 queryWithIdsConj, answer.queryPattern()
#         );
#
#         String expectedRule = explanationEntry.get("rule");
#         boolean hasExplanation = answer.hasExplanation();
#
#         if (expectedRule.equals("lookup")) {
#
#             assertFalse(String.format("Explanation entry %d is declared as a lookup, but an explanation was found", entryId), hasExplanation);
#
#             String[] expectedChildren = {"-"};
#             assertArrayEquals(String.format("Explanation entry %d is declared as a lookup, and so it should have no children, indicated as \"-\", but got children %s instead", entryId, Arrays.toString(children)), expectedChildren, children);
#         } else {
#
#             Explanation explanation = answer.explanation();
#             List<ConceptMap> explAnswers = explanation.getAnswers();
#
#             assertEquals(String.format("Explanation entry %d should have as many children as it has answers. Instead, %d children were declared, and %d answers were found.", entryId, children.length, explAnswers.size()), children.length, explAnswers.size());
#
#             if (expectedRule.equals("join")) {
#                 assertNull(String.format("Explanation entry %d is declared as a join, and should not have a rule attached, but one was found", entryId), explanation.getRule());
#             } else {
#                 // rule
#                 Rule.Remote rule = explanation.getRule();
#                 String ruleLabel = rule.getLabel();
#                 assertEquals(String.format("Incorrect rule label for explanation entry %d with rule %s.\nExpected: %s\nActual: %s", entryId, ruleLabel, expectedRule, ruleLabel), expectedRule, ruleLabel);
#
#                 Map<String, String> expectedRuleDefinition = rules.get(expectedRule);
#                 String when = Objects.requireNonNull(rule.getWhen()).toString();
#                 assertEquals(String.format("Incorrect rule body (when) for explanation entry %d with rule %s.\nExpected: %s\nActual: %s", entryId, ruleLabel, expectedRuleDefinition.get("when"), when), expectedRuleDefinition.get("when"), when);
#
#                 String then = Objects.requireNonNull(rule.getThen()).toString();
#                 assertEquals(String.format("Incorrect rule head (then) for explanation entry %d with rule %s.\nExpected: %s\nActual: %s", entryId, ruleLabel, expectedRuleDefinition.get("then"), then), expectedRuleDefinition.get("then"), then);
#             }
#             for (String child : children) {
#                 // Recurse
#                 checkExplanationEntry(explAnswers, explanationTree, Integer.valueOf(child));
#             }
#         }
#
#     @Then("each answer satisfies")
#     public void each_answer_satisfies(String templatedQuery) {
#         for (ConceptMap answer : answers) {
#             String query = applyQueryTemplate(templatedQuery, answer);
#             GraqlMatch graqlQuery = Graql.parseQuery(query).asMatch();
#             long answerSize = tx().query().match(graqlQuery).count();
#             assertEquals(1, answerSize);
#         }
#     }
#
#     private String applyQueryTemplate(String template, ConceptMap templateFiller) {
#         // find shortest matching strings between <>
#         Pattern pattern = Pattern.compile("<.+?>");
#         Matcher matcher = pattern.matcher(template);
#
#         StringBuilder builder = new StringBuilder();
#         int i = 0;
#         while (matcher.find()) {
#             String matched = matcher.group(0);
#             String requiredVariable = variableFromTemplatePlaceholder(matched.substring(1, matched.length() - 1));
#
#             builder.append(template, i, matcher.start());
#             if (templateFiller.map().containsKey(requiredVariable)) {
#
#                 Concept concept = templateFiller.get(requiredVariable);
#                 if (!concept.isThing())
#                     throw new ScenarioDefinitionException("Cannot apply IID templating to Type concepts");
#                 String conceptId = concept.asThing().getIID();
#                 builder.append(conceptId);
#
#             } else {
#                 throw new ScenarioDefinitionException(String.format("No IID available for template placeholder: %s.", matched));
#             }
#             i = matcher.end();
#         }
#         builder.append(template.substring(i));
#         return builder.toString();
#     }
#
#     private String variableFromTemplatePlaceholder(String placeholder) {
#         if (placeholder.endsWith(".iid")) {
#             String stripped = placeholder.replace(".iid", "");
#             String withoutPrefix = stripped.replace("answer.", "");
#             return withoutPrefix;
#         } else {
#             throw new ScenarioDefinitionException("Cannot replace template not based on ID.");
#         }
#     }
#
#     private static class ScenarioDefinitionException extends RuntimeException {
#         ScenarioDefinitionException(String message) {
#             super(message);
#         }
#     }
#
#     private interface UniquenessCheck {
#         boolean check(Concept concept);
#     }
#
#     public static class LabelUniquenessCheck implements UniquenessCheck {
#
#         private final Label label;
#
#         LabelUniquenessCheck(String scopedLabel) {
#             String[] tokens = scopedLabel.split(":");
#             this.label = tokens.length > 1 ? Label.of(tokens[0], tokens[1]) : Label.of(tokens[0]);
#         }
#
#         @Override
#         public boolean check(Concept concept) {
#             if (concept.isType()) return label.equals(concept.asType().getLabel());
#             throw new ScenarioDefinitionException("Concept was checked for label uniqueness, but it is not a Type.");
#         }
#     }
#
#     public static class AttributeUniquenessCheck {
#
#         protected final Label type;
#         protected final String value;
#
#         AttributeUniquenessCheck(String typeAndValue) {
#             String[] s = typeAndValue.split(":");
#             assertEquals(
#                     String.format("A check for attribute uniqueness should be given in the format \"type:value\", but received %s.", typeAndValue),
#                     2, s.length
#             );
#             type = Label.of(s[0]);
#             value = s[1];
#         }
#     }
#
#     public static class ValueUniquenessCheck extends AttributeUniquenessCheck implements UniquenessCheck {
#         ValueUniquenessCheck(String typeAndValue) {
#             super(typeAndValue);
#         }
#
#         public boolean check(Concept concept) {
#             if (!concept.isAttribute()) {
#                 return false;
#             }
#
#             Attribute<?> attribute = concept.asAttribute();
#             AttributeType attributeType = attribute.getType();
#
#             if (!type.equals(attributeType.getLabel())) {
#                 return false;
#             }
#
#             switch (attributeType.getValueType()) {
#                 case BOOLEAN:
#                     return value.equals(attribute.asBoolean().getValue().toString());
#                 case LONG:
#                     return value.equals(attribute.asLong().getValue().toString());
#                 case DOUBLE:
#                     return value.equals(attribute.asDouble().getValue().toString());
#                 case STRING:
#                     return value.equals(attribute.asString().getValue());
#                 case DATETIME:
#                     LocalDateTime dateTime;
#                     try {
#                         dateTime = LocalDateTime.parse(value);
#                     } catch (DateTimeParseException e) {
#                         dateTime = LocalDate.parse(value).atStartOfDay();
#                     }
#                     return dateTime.equals(attribute.asDateTime().getValue());
#                 case OBJECT:
#                 default:
#                     throw new ScenarioDefinitionException("Unrecognised value type " + attributeType.getValueType());
#             }
#         }
#     }
#
#     public static class KeyUniquenessCheck extends AttributeUniquenessCheck implements UniquenessCheck {
#         KeyUniquenessCheck(String typeAndValue) {
#             super(typeAndValue);
#         }
#
#         @Override
#         public boolean check(Concept concept) {
#             if (!concept.isThing()) { return false; }
#
#             Set<Attribute<?>> keys = concept.asThing().asRemote(tx()).getHas(true).collect(Collectors.toSet());
#             HashMap<Label, String> keyMap = new HashMap<>();
#
#             for (Attribute<?> key : keys) {
#                 String keyValue;
#                 switch (key.getType().getValueType()) {
#                     case BOOLEAN:
#                         keyValue = key.asBoolean().getValue().toString();
#                         break;
#                     case LONG:
#                         keyValue = key.asLong().getValue().toString();
#                         break;
#                     case DOUBLE:
#                         keyValue = key.asDouble().getValue().toString();
#                         break;
#                     case STRING:
#                         keyValue = key.asString().getValue();
#                         break;
#                     case DATETIME:
#                         keyValue = key.asDateTime().getValue().toString();
#                         break;
#                     case OBJECT:
#                     default:
#                         throw new ScenarioDefinitionException("Unrecognised value type " + key.getType().getValueType());
#                 }
#
#                 keyMap.put(key.getType().getLabel(), keyValue);
#             }
#             return value.equals(keyMap.get(type));
#         }
#     }
# }


text:
# SAM Record
# ==========

mutable struct Record
    # Data and filled range.
    data::Vector{UInt8}
    filled::UnitRange{Int} # Note: Specifies the data in use.

    # Mandatory fields.
    qname::UnitRange{Int}
    flag::UnitRange{Int}
    rname::UnitRange{Int}
    pos::UnitRange{Int}
    mapq::UnitRange{Int}
    cigar::UnitRange{Int}
    rnext::UnitRange{Int}
    pnext::UnitRange{Int}
    tlen::UnitRange{Int}
    seq::UnitRange{Int}
    qual::UnitRange{Int}

    # Auxiliary fields.
    fields::Vector{UnitRange{Int}}
end

"""
    SAM.Record()

Create an unfilled SAM record.
"""
function Record()
    return Record(
        UInt8[], 1:0,
        # qname-mapq
        1:0, 1:0, 1:0, 1:0, 1:0,
        # cigar-seq
        1:0, 1:0, 1:0, 1:0, 1:0,
        # qual and fields
        1:0, UnitRange{Int}[])
end

"""
    SAM.Record(data::Vector{UInt8})

Create a SAM record from `data`.
This function verifies the format and indexes fields for accessors.
Note that the ownership of `data` is transferred to a new record object.
"""
function Record(data::Vector{UInt8})
    return convert(Record, data)
end

function Base.convert(::Type{Record}, data::Vector{UInt8})
    record = Record(
        data, 1:0,
        # qname-mapq
        1:0, 1:0, 1:0, 1:0, 1:0,
        # cigar-seq
        1:0, 1:0, 1:0, 1:0, 1:0,
        # qual and fields
        1:0, UnitRange{Int}[])
    index!(record)
    return record
end

"""
    SAM.Record(str::AbstractString)

Create a SAM record from `str`.
This function verifies the format and indexes fields for accessors.
"""
function Record(str::AbstractString)
    return convert(Record, str)
end

function Base.convert(::Type{Record}, str::AbstractString)
    return Record(collect(UInt8, str))
end

function Base.:(==)(a::Record, b::Record)
    return a.filled == b.filled &&
        a.qname     == b.qname &&
        a.flag      == b.flag &&
        a.rname     == b.rname &&
        a.pos       == b.pos &&
        a.mapq      == b.mapq &&
        a.cigar     == b.cigar &&
        a.rnext     == b.rnext &&
        a.pnext     == b.pnext &&
        a.tlen      == b.tlen &&
        a.seq       == b.seq &&
        a.qual      == b.qual &&
        a.fields    == b.fields &&
        a.data[a.filled] == b.data[b.filled]
end

function Base.show(io::IO, record::Record)
    print(io, summary(record), ':')
    if isfilled(record)
        println(io)
        println(io, "    template name: ", hastempname(record) ? tempname(record) : "<missing>")
        println(io, "             flag: ", hasflag(record) ? flag(record) : "<missing>")
        println(io, "        reference: ", hasrefname(record) ? refname(record) : "<missing>")
        println(io, "         position: ", hasposition(record) ? position(record) : "<missing>")
        println(io, "  mapping quality: ", hasmappingquality(record) ? mappingquality(record) : "<missing>")
        println(io, "            CIGAR: ", hascigar(record) ? cigar(record) : "<missing>")
        println(io, "   next reference: ", hasnextrefname(record) ? nextrefname(record) : "<missing>")
        println(io, "    next position: ", hasnextposition(record) ? nextposition(record) : "<missing>")
        println(io, "  template length: ", hastemplength(record) ? templength(record) : "<missing>")
        println(io, "         sequence: ", hassequence(record) ? sequence(String, record) : "<missing>")
        println(io, "     base quality: ", hasquality(record) ? quality(String, record) : "<missing>")
          print(io, "   auxiliary data:")
        for field in record.fields
            print(io, ' ', String(record.data[field]))
        end
    else
        print(io, " <not filled>")
    end
end

function Base.print(io::IO, record::Record)
    write(io, record)
    return nothing
end

function Base.write(io::IO, record::Record)
    checkfilled(record)
    return unsafe_write(io, pointer(record.data, first(record.filled)), length(record.filled))
end

function Base.copy(record::Record)
    return Record(
        copy(record.data),
        record.filled,
        record.qname,
        record.flag,
        record.rname,
        record.pos,
        record.mapq,
        record.cigar,
        record.rnext,
        record.pnext,
        record.tlen,
        record.seq,
        record.qual,
        copy(record.fields))
end

# Accessor Functions
# ------------------

"""
    flag(record::Record)::UInt16

Get the bitwise flag of `record`.
"""
function flag(record::Record)::UInt16
    checkfilled(record)
    return unsafe_parse_decimal(UInt16, record.data, record.flag)
end

function hasflag(record::Record)
    return isfilled(record)
end

"""
    ismapped(record::Record)::Bool

Test if `record` is mapped.
"""
function ismapped(record::Record)::Bool
    return isfilled(record) && (flag(record) & FLAG_UNMAP == 0)
end

"""
    isprimary(record::Record)::Bool

Test if `record` is a primary line of the read.

This is equivalent to `flag(record) & 0x900 == 0`.
"""
function isprimary(record::Record)::Bool
    return flag(record) & 0x900 == 0
end

"""
    refname(record::Record)::String

Get the reference sequence name of `record`.
"""
function refname(record::Record)
    checkfilled(record)
    if ismissing(record, record.rname)
        missingerror(:refname)
    end
    return String(record.data[record.rname])
end

function hasrefname(record::Record)
    return isfilled(record) && !ismissing(record, record.rname)
end

"""
    position(record::Record)::Int

Get the 1-based leftmost mapping position of `record`.
"""
function position(record::Record)::Int
    checkfilled(record)
    pos = unsafe_parse_decimal(Int, record.data, record.pos)
    # if pos == 0
    #     missingerror(:position)
    # end
    return pos
end

function hasposition(record::Record)
    return isfilled(record) && (length(record.pos) != 1 || record.data[first(record.pos)] != UInt8('0'))
end

"""
    rightposition(record::Record)::Int

Get the 1-based rightmost mapping position of `record`.
"""
function rightposition(record::Record)
    return position(record) + alignlength(record) - 1
end

function hasrightposition(record::Record)
    return hasposition(record) && hasalignment(record)
end

"""
    isnextmapped(record::Record)::Bool

Test if the mate/next read of `record` is mapped.
"""
function isnextmapped(record::Record)::Bool
    return isfilled(record) && (flag(record) & FLAG_MUNMAP == 0)
end

"""
    nextrefname(record::Record)::String

Get the reference name of the mate/next read of `record`.
"""
function nextrefname(record::Record)::String
    checkfilled(record)
    if ismissing(record, record.rnext)
        missingerror(:nextrefname)
    end
    return String(record.data[record.rnext])
end

function hasnextrefname(record::Record)
    return isfilled(record) && !ismissing(record, record.rnext)
end

"""
    nextposition(record::Record)::Int

Get the position of the mate/next read of `record`.
"""
function nextposition(record::Record)::Int
    checkfilled(record)
    pos = unsafe_parse_decimal(Int, record.data, record.pnext)
    # if pos == 0
    #     missingerror(:nextposition)
    # end
    return pos
end

function hasnextposition(record::Record)
    return isfilled(record) && (length(record.pnext) != 1 || record.data[first(record.pnext)] != UInt8('0'))
end

"""
    mappingquality(record::Record)::UInt8

Get the mapping quality of `record`.
"""
function mappingquality(record::Record)::UInt8
    checkfilled(record)
    qual = unsafe_parse_decimal(UInt8, record.data, record.mapq)
    if qual == 0xff
        missingerror(:mappingquality)
    end
    return qual
end

function hasmappingquality(record::Record)
    return isfilled(record) && unsafe_parse_decimal(UInt8, record.data, record.mapq) != 0xff
end

"""
    cigar(record::Record)::String

Get the CIGAR string of `record`.
"""
function cigar(record::Record)::String
    checkfilled(record)
    if ismissing(record, record.cigar)
        # missingerror(:cigar)
        return ""
    end
    return String(record.data[record.cigar])
end

function hascigar(record::Record)
    return isfilled(record) && !ismissing(record, record.cigar)
end

"""
    alignment(record::Record)::BioAlignments.Alignment

Get the alignment of `record`.
"""
function alignment(record::Record)::BioAlignments.Alignment
    if ismapped(record)
        return BioAlignments.Alignment(cigar(record), 1, position(record))
    end

    return BioAlignments.Alignment(BioAlignments.AlignmentAnchor[])

end

function hasalignment(record::Record)
    return isfilled(record) && hascigar(record)
end

"""
    alignlength(record::Record)::Int

Get the alignment length of `record`.
"""
function alignlength(record::Record)::Int
    if length(record.cigar) == 1 && record.data[first(record.cigar)] == UInt8('*')
        return 0
    end
    ret::Int = 0
    len = 0  # operation length
    for i in record.cigar
        c = record.data[i]
        if in(c, UInt8('0'):UInt8('9'))
            len = len * 10 + (c - UInt8('0'))
            continue
        end
        op = convert(BioAlignments.Operation, Char(c))
        if BioAlignments.ismatchop(op) || BioAlignments.isdeleteop(op) #Note: reference consuming ops ('M', 'D', 'N', '=', 'X').
            ret += len
        end
        len = 0
    end
    return ret
end

"""
    tempname(record::Record)::String

Get the query template name of `record`.
"""
function tempname(record::Record)::String
    checkfilled(record)
    if ismissing(record, record.qname)
        missingerror(:tempname)
    end
    return String(record.data[record.qname])
end

function hastempname(record::Record)
    return isfilled(record) && !ismissing(record, record.qname)
end

"""
    templength(record::Record)::Int

Get the template length of `record`.
"""
function templength(record::Record)::Int
    checkfilled(record)
    len = unsafe_parse_decimal(Int, record.data, record.tlen)
    # if len == 0
    #     missingerror(:tlen)
    # end
    return len
end

function hastemplength(record::Record)
    return isfilled(record) && (length(record.tlen) != 1 || record.data[first(record.tlen)] != UInt8('0'))
end

"""
    sequence(record::Record)::BioSequences.LongDNASeq

Get the segment sequence of `record`.
"""
function sequence(record::Record)
    checkfilled(record)
    if ismissing(record, record.seq)
        # missingerror(:sequence)
        return nothing
    end
    seqlen = length(record.seq)
    ret = BioSequences.LongDNASeq(seqlen)
    copyto!(ret, 1, record.data, first(record.seq), seqlen)
    return ret
end

function hassequence(record::Record)
    return isfilled(record) && !ismissing(record, record.seq)
end
"""
    sequence(::Type{String}, record::Record)::String

Get the segment sequence of `record` as `String`.
"""
function sequence(::Type{String}, record::Record)::String
    checkfilled(record)
    return String(record.data[record.seq])
end

"""
    seqlength(record::Record)::Int

Get the sequence length of `record`.
"""
function seqlength(record::Record)::Int
    checkfilled(record)
    if ismissing(record, record.seq)
        missingerror(:seq)
    end
    return length(record.seq)
end

function hasseqlength(record::Record)
    return isfilled(record) && !ismissing(record, record.seq)
end

"""
    quality(record::Record)::Vector{UInt8}

Get the Phred-scaled base quality of `record`.
"""
function quality(record::Record)::Vector{UInt8}
    checkfilled(record)
    if ismissing(record, record.qual)
        missingerror(:quality)
    end
    qual = record.data[record.qual]
    for i in 1:lastindex(qual)
        @inbounds qual[i] -= 33
    end
    return qual
end

function hasquality(record::Record)
    return isfilled(record) && !ismissing(record, record.qual)
end

"""
    quality(::Type{String}, record::Record)::String

Get the ASCII-encoded base quality of `record`.
"""
function quality(::Type{String}, record::Record)::String
    checkfilled(record)
    return String(record.data[record.qual])
end

"""
    auxdata(record::Record)::Dict{String,Any}

Get the auxiliary data (optional fields) of `record`.
"""
function auxdata(record::Record)::Dict{String,Any}
    checkfilled(record)
    return Dict(k => record[k] for k in keys(record))
end

function Base.haskey(record::Record, tag::AbstractString)
    return findauxtag(record, tag) > 0
end

function Base.getindex(record::Record, tag::AbstractString)
    i = findauxtag(record, tag)
    if i == 0
        throw(KeyError(tag))
    end
    field = record.fields[i]
    # data type
    typ = record.data[first(field)+3]
    lo = first(field) + 5
    if i == lastindex(record.fields)
        hi = last(field)
    else
        hi = first(record.fields[i+1]) - 2
    end

    if typ == UInt8('A')
        @assert lo == hi
        return Char(record.data[lo])
    end
    if typ == UInt8('i')
        return unsafe_parse_decimal(Int, record.data, lo:hi)
    end
    if typ == UInt8('f')
        # TODO: Call a C function directly for speed?
        return parse(Float32, SubString(record.data[lo:hi]))
    end
    if typ == UInt8('Z')
        return String(record.data[lo:hi])
    end
    if typ == UInt8('H')
        return parse_hexarray(record.data, lo:hi)
    end
    if typ == UInt8('B')
        return parse_typedarray(record.data, lo:hi)
    end

    throw(ArgumentError("type code '$(Char(typ))' is not defined"))

end

function Base.keys(record::Record)
    checkfilled(record)
    return [String(record.data[first(f):first(f)+1]) for f in record.fields]
end

function Base.values(record::Record)
    return [record[k] for k in keys(record)]
end

# Bio Methods
# -----------

function BioGenerics.isfilled(record::Record)
    return !isempty(record.filled)
end

function BioGenerics.seqname(record::Record)
    return tempname(record)
end

function BioGenerics.hasseqname(record::Record)
    return hastempname(record)
end

function BioGenerics.sequence(record::Record)
    return sequence(record)
end

function BioGenerics.hassequence(record::Record)
    return hassequence(record)
end

function BioGenerics.rightposition(record::Record)
    return rightposition(record)
end

function BioGenerics.hasrightposition(record::Record)
    return hasrightposition(record)
end

function BioGenerics.leftposition(record::Record)
    return position(record)
end

function BioGenerics.hasleftposition(record::Record)
    return hasposition(record)
end

# Helper Functions
# ----------------

function Base.empty!(record::Record)
    record.filled = 1:0
    record.qname = 1:0
    record.flag = 1:0
    record.rname = 1:0
    record.pos = 1:0
    record.mapq = 1:0
    record.cigar = 1:0
    record.rnext = 1:0
    record.pnext = 1:0
    record.tlen = 1:0
    record.seq = 1:0
    record.qual = 1:0
    empty!(record.fields)
    return record
end

function initialize!(record::Record) #TODO: deprecate.
    return empty!(record)
end

function checkfilled(record::Record)
    if !isfilled(record)
        throw(ArgumentError("unfilled SAM record"))
    end
end

function findauxtag(record::Record, tag::AbstractString)
    checkfilled(record)
    if sizeof(tag) != 2
        return 0
    end
    t1, t2 = UInt8(tag[1]), UInt8(tag[2])
    for (i, field) in enumerate(record.fields)
        p = first(field)
        if record.data[p] == t1 && record.data[p+1] == t2
            return i
        end
    end
    return 0
end

function parse_hexarray(data::Vector{UInt8}, range::UnitRange{Int})
    @assert iseven(length(range))
    ret = Vector{UInt8}(length(range) >> 1)
    byte2hex(b) = b ∈ 0x30:0x39 ? (b - 0x30) : b ∈ 0x41:0x46 ? (b - 0x41 + 0x0A) : error("not in [0-9A-F]")
    j = 1
    for i in first(range):2:last(range)-1
        ret[j] = (byte2hex(data[range[i]]) << 4) | byte2hex(data[range[i+1]])
        j += 1
    end
    return ret
end

function parse_typedarray(data::Vector{UInt8}, range::UnitRange{Int})
    # format: [cCsSiIf](,[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)+
    t = data[first(range)]
    xs = split(String(data[first(range)+2:last(range)]))
    if t == UInt8('c')
        return [parse(Int8, x) for x in xs]
    end
    if t == UInt8('C')
        return [parse(UInt8, x) for x in xs]
    end
    if t == UInt8('s')
        return [parse(Int16, x) for x in xs]
    end
    if t == UInt8('S')
        return [parse(UInt16, x) for x in xs]
    end
    if t == UInt8('i')
        return [parse(Int32, x) for x in xs]
    end
    if t == UInt8('I')
        return [parse(UInt32, x) for x in xs]
    end
    if t == UInt8('f')
        return [parse(Float32, x) for x in xs]
    end

    throw(ArgumentError("type code '$(Char(t))' is not defined"))

end

function ismissing(record::Record, range::UnitRange{Int})
    return length(range) == 1 && record.data[first(range)] == UInt8('*')
end

ext: jl


text:
using Test
using Bridge
X = SamplePath(0.0:1.0:10.0,Bridge.cumsum0(rand(10)))

#@test_throws SamplePath(0.0:1.0:10.0,Bridge.cumsum0(rand(10)))

@test length(X) == length(X.tt)

@test collect(t for (t,x) in X) == X.tt
@test collect(x for (t,x) in X) == X.yy

@test Base.iteratorsize(SamplePath) == Base.HasLength()


text:
module Tests

include("utilities.jl")

const linear_tests = Dict{String, Function}()

import MultivariateBases
const MB = MultivariateBases

include("zero_polynomial.jl")
include("zero_polynomial_in_fixed_variables_set.jl")
include("zero_polynomial_in_algebraic_set.jl")
include("plus_minus.jl")

@test_suite linear

end


text:
module Piecewise

    using Measurements

    include("common.jl")
    include(joinpath("staticpolynomials","StaticPolynomials.jl"))
    include(joinpath("staticpolynomials","parsing.jl"))

    include("piecewisepolynomials.jl")
    include("orderedpiecewise.jl")
    include("macros.jl")

 end


text:
using Pkg, Revise, FourierGPE
gr(titlefontsize=12,size=(500,300),transpose=true,colorbar=false)

## set simulation parameters
L = (25.0,25.0)
N = (256,256)
μ = 25.0

## Initialize simulation
sim = Sim(L,N)
@pack! sim = μ
@unpack_Sim sim

## Two ways to set potential:

# Time dependent potential function (here trivial t dep)
import FourierGPE.V
V(x,y,t) = 0.5*(x^2 + y^2)
ψ0(x,y,μ,g) = sqrt(μ/g)*sqrt(max(1.0-V(x,y,0.0)/μ,0.0)+im*0.0)

# As a static Potential
# x,y = X
# Vs(x,y) = 0.5*(x^2 + y^2)
# ψ0(x,y,μ,g) = sqrt(μ/g)*sqrt(max(1.0-Vs(x,y)/μ,0.0)+im*0.0)
# sim = Sim(sim,V0=Vs.(x,y'))

## make initial state
x,y = X
ψi = ψ0.(x,y',μ,g)
ϕi = kspace(ψi,sim)
@pack! sim = ϕi,μ

## Evolve in k space
@time sol = runsim(sim)

## ground state
ϕg = sol[end]
ψg = xspace(ϕg,sim)
showpsi(x,y,ψg)

## set simulation parameters
γ = 0.3
t = LinRange(ti,tf,Nt)
ϕi = kspace(ψg,sim)
# reltol = 1e-7
# alg = DP5()
@pack! sim = tf,t,γ,ϕi #,reltol,alg
@unpack_Sim sim

import FourierGPE: nlin!, Lgp!
sim

function Lz!(dϕ,ϕ,sim::Sim{2},t)
    #TODO
    @unpack X,K,T,params = sim; x,y = X; kx, ky = K
    χ = T.psi.x[1]
    Ω = params.κ

    # x*ky
    χ .= ϕ
    @. χ *= ky'
    xspace!(χ,sim)
    @. χ *= x
    kspace!(χ,sim)
    @. dϕ += Ω*χ

    # -y*kx
    χ .= ϕ
    @. χ *= kx
    xspace!(χ,sim)
    @. χ *= y'
    kspace!(χ,sim)
    @. dϕ -= Ω*χ
    return nothing
end

function Lz(ϕ,sim::Sim{2})
    @unpack X,K,params = sim; x,y = X; kx, ky = K
    Ω = params.κ
    # x*ky
    psi1 = ϕ.*ky'
    xspace!(psi1,sim)
    @. psi1 *= Ω*x
    kspace!(psi1,sim)

    # y*kx
    psi2 = ϕ.*kx
    xspace!(psi2,sim)
    @. psi2 *= Ω*y'
    kspace!(psi2,sim)
    return psi1 .- psi2
end

function nlin!(dϕ,ϕ,sim::Sim{2},t)
    @unpack g,X,V0 = sim; x,y = X
    y = y'
    dϕ .= ϕ
    xspace!(dϕ,sim)
    @. dϕ *= V0 + V(x,y,t) + g*abs2(dϕ)
    kspace!(dϕ,sim)
    return nothing
end

function Lgp!(dϕ,ϕ,sim,t)
    @unpack γ,μ,espec = sim
    nlin!(dϕ,ϕ,sim,t)
    dϕ .+= Lz(ϕ,sim)
    @. dϕ = -im*(1.0 - im*γ)*(dϕ + (espec - μ)*ϕ)
    return nothing
end

params = Params(κ = 0.9)
@pack! sim = params
sim
ti = 0.; tf = 5.
t = LinRange(ti,tf,Nt)
ϕi = ϕi + randn(N)+im*randn(N)
# ==== Evolve in k space
solv = runsim(sim)

# ===================================
ϕf = solv[end]
ψf = xspace(ϕf,sim)
showpsi(x,y,ψf)

anim = @animate for i in eachindex(t)
    ψ = xspace(solv[i],sim)
    showpsi(x,y,ψ)
end

gif(anim,"./examples/vortexlattice.gif",fps=25)


text:
# This file is part of IntegerTriangles.jl.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module TrianglesBase

using Nemo
import Base.sum

export AbstractSequence, AbstractTriangle
export ℤInt, ℤSeq, ℤTri, ℚInt, ℚSeq, ℚTri
export ℤPolySeq, ℚPolySeq, ℤPolyRing, ℚPolyRing
export ZInt, ZSeq, ZTri, QInt, QSeq, QTri
export ZPolySeq, QPolySeq, ZPolyRing, QPolyRing
export Polynomial, Evaluate, PolyArray
export PolyTriangle, PolyTri, DiagTri, DiagonalTriangle
export OrthoPoly, DelehamΔ, RecTriangle
export EgfExpansionCoeff, EgfExpansionPoly
export InverseTriangle, Inverse, Reverse, InvRev, RevInv
export PolynomialFunction, ReversePolynomial, PolyValue
export PolyVal2, PolyVal3, PolyValn, Coefficients
export RiordanSquare, Println, Print, Telescope, LinMap, Trans
export EvenSum, OddSum, AltSum, DiagSum, Central, Middle
export LeftSide, RightSide, PosHalf, NegHalf, Flat, Factorial
export Binomial, BinomialTransform, BinConv, BinomialTriangle
export IBinomial, IBinomialTransform, IBinConv, IBinomialTriangle
export TransUnos, TransAlts, TransSqrs, TransNat0, TransNat1
export Convolution, ConvolutionTriangle, ConvolutionTransformation
export ConvSum, ConvTri, ConvTrans
export TRAITS

"""
Basic definitions used in the library IntegerTriangles.
We consider the following traits of an triangle:

```
Flat, Reverse, Inverse, RevInv, InvRev, DiagTri,

PolyTri, Sum, EvenSum, OddSum, AltSum, DiagSum,

Middle, Central, LeftSide, RightSide, PosHalf,

NegHalf, PolyVal2, PolyVal3, BinConv, IBinConv,

TransSqrs, TransNat0, TransNat1
```
"""
const ModuleTrianglesBase = ""

"""
Supertype for sequences (or sequence-like types).
"""
abstract type AbstractSequence end

"""
Supertype for triangles (or triangles-like types).
"""
abstract type AbstractTriangle end

const Seq{T} = Array{T,1}

"""
ℤInt is an alias for the type Nemo.fmpz.
"""
const ℤInt  = Nemo.fmpz

"""
ℚInt is an alias for the type Nemo.fmpq.
"""
const ℚInt  = Nemo.fmpq

"""
ℤSeq is an alias for the type Array{Nemo.fmpz, 1}.
"""
const ℤSeq = Seq{ℤInt}

"""
ℚSeq is an alias for the type Array{Nemo.fmpq, 1}.
"""
const ℚSeq = Seq{ℚInt}

"""
ℤTri is an alias for the type Array{Array{Nemo.fmpz, 1}, 1}.
"""
const ℤTri = Seq{ℤSeq}

"""
ℚTri is an alias for the type Array{Array{Nemo.fmpq, 1}, 1}.
"""
const ℚTri = Seq{ℚSeq}

# -------------

# ZZ(n)    # (constructor for fmpz, defined in Nemo)
# QQ(n, k) # (constructor for fmpq, defined in Nemo)

"""
Constructor for an ℤSeq of length len. If a second parameter f
is given the sequence is constructed as [ZZ(f(n)) for n in 0:len-1]
"""
ZSeq(len::Int) = ℤSeq(undef, len)

ZSeq(len::Int, f::Function) = [ZZ(f(n)) for n in 0:len-1]

"""
Constructor for an ℚSeq of length len. If a second parameter f
is given the sequence is constructed as [QQ(f(n)) for n in 0:len-1]
"""
QSeq(len::Int) = ℚSeq(undef, len)

QSeq(len::Int, f::Function) = [QQ(f(n)) for n in 0:len-1]

"""
Constructor for ZTri.
"""
function ZTri(dim::Int; reg=false)
    reg ? ℤTri(ZSeq.(1:dim)) : ℤTri(undef, dim)
end

ZTri(dim, f::Function) = f.(0:dim-1)
ZTri(T::AbstractTriangle) = [row for row ∈ T]

"""
Constructor for QTri.
"""
function QTri(dim::Int; reg=false)
    reg ? QSeq.(1:dim) : ℚTri(undef, dim)
end

# ---------------

"""
ℤPoly is an alias for the type Nemo.fmpz_poly.
"""
const ℤPoly = Nemo.fmpz_poly

"""
ℚPoly is an alias for the type Nemo.fmpq_poly.
"""
const ℚPoly = Nemo.fmpq_poly

"""
ℤPolySeq is an alias for the type Array{Nemo.fmpz_poly, 1}.
"""
const ℤPolySeq = Seq{ℤPoly}

"""
ℚPolySeq is an alias for the type Array{Nemo.fmpq_poly, 1}.
"""
const ℚPolySeq = Seq{ℚPoly}

"""
ℤPolyTri is an alias for the type Array{Array{Nemo.fmpz_poly, 1}, 1}.
"""
const ℤPolyTri = Seq{ℤPolySeq}

"""
ℚPolyTri is an alias for the type Array{Array{Nemo.fmpq_poly, 1}, 1}.
"""
const ℚPolyTri = Seq{ℚPolySeq}

# ---------------

"""
(SIGNATURES)
"""
ZPolyRing(x) = PolynomialRing(ZZ, x)

"""
(SIGNATURES)
"""
QPolyRing(x) = PolynomialRing(QQ, x)

const WARNING_ON_NOTINVERTIBLE = false

"""
Return the diagonal triangle T(n - k, k) where k in 0:n/2.
```
julia> Println.(DiagonalTriangle(MotzkinTriangle(8)))
[1]
[1]
[2, 1]
[4, 2]
[9, 5, 1]
[21, 12, 3]
[51, 30, 9, 1]
[127, 76, 25, 4]  # A106489
```
"""
function DiagonalTriangle(T::ℤTri)
    dim = length(T)
    U = ZTri(dim)
    for n = 1:dim
        R = ZSeq(div(n + 1, 2))
        for k = 0:div(n - 1, 2)
            R[k + 1] = T[n - k][k + 1]
        end
        U[n] = R
    end
    U
end

DiagTri(T::ℤTri) = DiagonalTriangle(T)

# --------------------------------------------------

function Convolution(S::ℤSeq, T::ℤSeq)
    L = length(S);
    if L != length(T)
        @warn "Incompatible lengths!"
        return []
    end
    [sum(S[k + 1] * T[n - k + 1] for k in 0:n) for n in 0:L-1]
end

Convolution(S::ℤSeq) = Convolution(S, S)
ConvSum(S::ℤSeq) = sum(Convolution(S))

"""
Return the convolution triangle of T.
```
julia> T = [ℤInt[1], ℤInt[2, 3], ℤInt[4, 5, 6], ℤInt[7, 8, 9, 10]]
Println.(ConvolutionTriangle(T))
[1]
[4, 12]
[16, 40, 73]
[49, 112, 190, 284]
```
"""
ConvolutionTriangle(T::ℤTri) = Convolution.(T)
ConvTri(T::ℤTri) = ConvolutionTriangle(T) # alias

#= --------------------------------------------------

"""
```
julia> Println.(ConvolutionTransformation(T))
[1, 16, 129, 635]
```
"""
ConvolutionTransformation(T::ℤTri) = sum.(Convolution.(T))
ConvTrans(T::ℤTri) = ConvolutionTransformation(T::ℤTri) # alias

function BinomialConvolution(S::ℤSeq, T::ℤSeq)
    L = length(S); L != length(T) && return []
    [sum(Binomial(n,k)*S[k + 1] * T[n - k + 1] for k in 0:n) for n in 0:L-1]
end

BinomialConvolution(S::ℤSeq) = BinomialConvolution(S, S)
BinConvSum(S::ℤSeq) = sum(BinomialConvolution(S))

"""
Return the binomial convolution triangle of T.
"""
BinomialConvolutionTriangle(T::ℤTri) = BinomialConvolution.(T)
BinConvTri(T::ℤTri) = BinomialConvolutionTriangle(T) # alias

BinomialConvolutionTransformation(T::ℤTri) = sum.(BinomialConvolution.(T))
BinConvTrans(T::ℤTri) = BinomialConvolutionTransformation(T::ℤTri) # alias

=#

#--------------------------------------------------

"""
The sum of a ℤTri is the sequence of the sum of the rows.
"""
Base.sum(T::ℤTri) = [sum(row) for row ∈ T]
Base.sum(T::ℚTri) = [sum(row) for row ∈ T]

"""
The EvenSum of a ℤSeq is the sum of the even indexed terms, indexing starts with 0.
```
julia> EvenSum([0, 1, 2, 3, 4, 5])
0 + 2 + 4 = 6
```
"""
EvenSum(A) = sum(A[1:2:end])

"""
The OddSum of a ℤSeq is the sum of the odd indexed terms, indexing starts with 0.
```
julia> OddSum([0, 1, 2, 3, 4, 5])
1 + 3 + 5 = 9
```
"""
OddSum(A) = sum(A[2:2:end])

"""
The AltSum of a ℤSeq is the alternating sum.
```
julia> AltSum([0, 1, 2, 3, 4, 5])
+ 0 - 1 + 2 - 3 + 4 - 5 = 6 - 9 = - 3
```
"""
AltSum(A) = EvenSum(A) - OddSum(A)

"""
The Middle of a ℤSeq A is the middle term, A[div(begin + end, 2)].
```
julia> Middle([0, 1, 2, 3, 4, 5])
3
```
"""
Middle(A) = A[div(begin + end, 2)]

"""
The LeftSide of a ℤSeq is the first term.
```
julia> LeftSide([0, 1, 2, 3, 4, 5])
0
```
"""
LeftSide(A) = A[begin]

"""
The RightSide of a ℤSeq is the last term.
```
julia> RightSide([0, 1, 2, 3, 4, 5])
5
```
"""
RightSide(A) = A[end]

"""
The EvenSum of a ℤTri is the sequence of the even sums of the rows.
```
EvenSum(T::ℤTri) = EvenSum.(T)
```
"""
EvenSum(T::ℤTri) = EvenSum.(T)

"""
The OddSum of a ℤTri is the sequence of the odd sums of the rows.
```
OddSum(T::ℤTri) = OddSum.(T)
```
"""
OddSum(T::ℤTri) = OddSum.(T)

"""
The AltSum of a ℤTri is the sequence of the alternating sums of the rows.
```
AltSum(T::ℤTri) = EvenSum(T) - OddSum(T)
```
"""
AltSum(T::ℤTri) = EvenSum(T) - OddSum(T)

"""
The DiagSum of a ℤTri is the sum of the diagonal triangle.
```
DiagSum(T::ℤTri) = sum(DiagonalTriangle(T))
```
"""
DiagSum(T::ℤTri) = sum(DiagonalTriangle(T))

"""
The Middle of a ℤTri is the sequence of the middle term in the rows.
```
Middle(T::ℤTri) = Middle.(T)
```
"""
Middle(T::ℤTri) = Middle.(T)

"""
The Central of a ℤTri is the sequence of the middle term
of the even indexed rows, indexing starts with 0.
```
Central(T::ℤTri) = Middle.(T[begin:2:end])
```
"""
Central(T::ℤTri) = Middle.(T[begin:2:end])

"""
The LeftSide of a ℤTri is the sequence of the first term in the rows.
```
LeftSide(T::ℤTri) = LeftSide.(T)
```
"""
LeftSide(T::ℤTri) = LeftSide.(T)

"""
The RightSide of a ℤTri is the sequence of the last term in the rows.
```
RightSide(T::ℤTri) = RightSide.(T)
```
"""
RightSide(T::ℤTri) = RightSide.(T)

"""
Return the ℤ-polynomial whose coefficients are the terms of the sequence.
```
[1, 2, 3] -> 1 + 2*x + 3*x^2
```
"""
function Polynomial(S::ℤSeq)
    R, x = ZPolyRing("x")
    sum(c * x^(k - 1) for (k, c) ∈ enumerate(S))
end

"""
Return the sequence of ℤ-polynomials whose coefficients are the terms of the triangle.
```
Polynomial(T::ℤTri) = Polynomial.(T)
```
"""
Polynomial(T::ℤTri) = Polynomial.(T)

"""
Return the ℤ-polynomial whose coefficients are the terms of the reversed sequence.
```
[1, 2, 3] -> x^2 + 2*x + 3
```
"""
ReversePolynomial(S::ℤSeq) = Polynomial(reverse(S))

"""
Return the sequence of ℤ-polynomials whose coefficients are the terms of the reversed triangle.
```
ReversePolynomial(T::ℤTri) = ReversePolynomial.(T)
```
"""
ReversePolynomial(T::ℤTri) = ReversePolynomial.(T)

"""
Return the polynomial function associated with the polynomial with coefficients
given by the sequence S. A polynomial function evaluates to Float64 values.
```
p = PolynomialFunction([1, 2, 3])
julia> [p(r) for r in 0:3]
4-element Vector{Float64}:
  1.0
  6.0
 17.0
 34.0
```
"""
function PolynomialFunction(S)
    y -> sum(Float64(c) * y^(k - 1)
            for (k, c) ∈ enumerate(S))
end

"""
Return the sequence of polynomial functions associated to the polynomials with coefficients
given by the triangle T.
```
PolynomialFunction(T::ℤTri) = PolynomialFunction.(T)
```
"""
PolynomialFunction(T::ℤTri) = PolynomialFunction.(T)

"""
Return the value of the ℤ-polynomial p evaluated at x.
```
julia> R, x = ZPolyRing("x")
p = 1 + 2*x + 3*x^2
Evaluate(p, 2)
17
```
"""
Evaluate(p::ℤPoly, x) = subst(p, x)
Evaluate(p::ℚPoly, x) = subst(p, x)

"""
Return the sequence of values to which the sequence of ℤ-polynomials P evaluate at x.
```
julia> R, x = ZPolyRing("x")
P = [sum(k * x^k for k in 0:n) for n in 0:9]
Evaluate(P, 3)  # A289399
```
"""
Evaluate(P::ℤPolySeq, x) = Evaluate.(P, x)
Evaluate(P::ℚPolySeq, x) = Evaluate.(P, x)

"""
Return the list of coefficients of the polynomial p (ordered by ascending powers).
```
julia> R, x = ZPolyRing("x")
p = 1 + 2*x + 3*x^2
Coefficients(p)
```
"""
Coefficients(p, n=0) = p == 0 ? ℤInt[0] : coeff.(p, 0:max(n,degree(p)))

"""
Return the sequence of list of coefficients of the polynomials P.
"""
Coefficients(P::AbstractArray) = Coefficients.(P)

import Nemo.numerator
numerator(T::ℚTri) = [numerator.(t) for t ∈ T]

"""
Return numerator(2^degree(p) * Evaluate(p, QQ(1, 2))).
```
julia> R, x = ZPolyRing("x")
p = 1 + 2*x + 3*x^2
PosHalf(p)
11
```
"""
PosHalf(p) = p == 0 ? ZZ(0) : numerator(2^degree(p) * Evaluate(p, QQ(1, 2)))

"""
Return Numerator((-2)^degree(p) * Evaluate(p, QQ(-1, 2)))
```
julia> R, x = ZPolyRing("x")
p = 1 + 2*x + 3*x^2
NegHalf(p)
3
```
"""
NegHalf(p) = p == 0 ? ZZ(0) : numerator((-2)^degree(p) * Evaluate(p, QQ(-1, 2)))

"""
Return the sequence generated by broadcasting PosHalf over the
rows of the triangle interpreted as a polynomials.
```
julia> T = [[ZZ(k) for k in 0:n] for n in 1:5]
PosHalf(Polynomial(T))  # A000295
5-element ℤSeq
 1
 4
 11
 26
 57
```
"""
PosHalf(T::ℤTri) = PosHalf.(Polynomial(T))

"""
Return the sequence generated by broadcasting NegHalf over the
rows of the triangle interpreted as a polynomials.
```
NegHalf(T::ℤTri) = NegHalf.(Polynomial(T))
```
"""
NegHalf(T::ℤTri) = NegHalf.(Polynomial(T))

"""
Return the sequence generated by broadcasting PosHalf over the
sequence of polynomials.
"""
PosHalf(P::ℤPolySeq) = PosHalf.(P)

"""
Return the sequence generated by broadcasting NegHalf over the
sequence of polynomials.
"""
NegHalf(P::ℤPolySeq) = NegHalf.(P)

# -------------------------------------------------------

LinMap(M::ℤTri, V::ℤSeq, n) = sum(M[n][k] * V[k] for k = 1:n)

"""
LinMap(M::ℤTri, V::ℤSeq) returns the 'matrix times vector' product
of M and V. Expands up to length(V) and we assume length(M) >= length(V).
"""
LinMap(M::ℤTri, V::ℤSeq) = (n -> LinMap(M, V, n)).(1:length(V))

LinMap(F::Function, V::ℤSeq, n) = sum(F(n - 1)[k] * V[k] for k = 1:n)

"""
LinMap(F::Function, V::ℤSeq) returns the 'matrix times vector' product
of a matrix generated by F and V.
```
julia>  L(n) = [ZZ(k) for k in 0:n]
[LinMap(L, L(n)) for n in 0:9]
0
1
5
14
30
... # A000330
```
"""
LinMap(F::Function, V::ℤSeq) = LinMap(F, V, length(V))

"""
Transform a ℤSeq V by the matrix/vector product by applying LinMap to (M, V).
Expands up to min(length(M), length(V)).
"""
Trans(M::ℤTri, V::ℤSeq) = (n -> LinMap(M, V, n)).(1:min(length(M), length(V)))

# -------------------------------------------------------

#"""  --- same as sum
#TransUnos(T) = Trans(T, [ZZ(1) for n = 0:length(T)])
#
#Standard trait of T.
#"""
#TransUnos(T) = Trans(T, [ZZ(1) for n = 0:length(T)])

"""
TransAlts(T) = Trans(T, [(-1)^n * ZZ(1) for n = 0:length(T)])

Standard trait of T.
"""
TransAlts(T) = Trans(T, [(-1)^n * ZZ(1) for n = 0:length(T)])

"""
TransSqrs(T) = Trans(T, [ZZ(n^2) for n = 0:length(T)])

Standard trait of T.
"""
TransSqrs(T) = Trans(T, [ZZ(n^2) for n = 0:length(T)])

"""
TransNat0(T) = Trans(T, [ZZ(n) for n = 0:length(T)])

Standard trait of T.
"""
TransNat0(T) = Trans(T, [ZZ(n) for n = 0:length(T)])

"""
TransNat1(T) = Trans(T, [ZZ(n) for n = 1:length(T)])

Standard trait of T.
"""
TransNat1(T) = Trans(T, [ZZ(n) for n = 1:length(T)])

"""
Return an iterator expanding the given sequence to a regular triangle.
```
julia> T = Telescope(ℤInt[1, 2, 3, 4, 5, 6])
collect(T)
 [1]
 [1, 2]
 [1, 2, 3]
 [1, 2, 3, 4]
 [1, 2, 3, 4, 5]
 [1, 2, 3, 4, 5, 6]
```
"""
Telescope(A::ℤSeq) = (A[1:k] for k = 1:size(A, 1))

"""
Return an iterator expanding the sequence generated by the function f to a regular triangle.
"""
Telescope(len::Int, f::Function) = Telescope([ZZ(f(k)) for k = 0:len - 1])

# --------------------------------------------------

"""
(SIGNATURES)
"""
Binomial(n, k) = Nemo.binomial(ZZ(n), ZZ(k))
Binomial(n) = [Binomial(n, k) for k = 0:n]
Binomial(A::ℤSeq) = LinMap(Binomial, A)

"""
(SIGNATURES)
"""
BinomialTriangle(dim) = [Binomial(n) for n = 0:dim - 1]
"""
(SIGNATURES)
"""
BinomialTransform(A::ℤSeq) = Binomial.(Telescope(A))
BinomialTransform(A::ℤTri) = Binomial.(A)

# Alias
BinConv(A::ℤTri) = BinomialTransform(A::ℤTri)

# --------------------------------------------------

IBinomial(n, k) = ZZ((-1)^(n-k))*Binomial(n, k)
IBinomial(n) = [IBinomial(n, k) for k = 0:n]
IBinomial(A::ℤSeq) = LinMap(IBinomial, A)

"""
(SIGNATURES)
"""
IBinomialTriangle(dim) = [IBinomial(n) for n = 0:dim - 1]
"""
(SIGNATURES)
"""
IBinomialTransform(A::ℤSeq) = IBinomial.(Telescope(A))
IBinomialTransform(A::ℤTri) = IBinomial.(A)

# Alias
IBinConv(A::ℤTri) = IBinomialTransform(A::ℤTri)

# --------------------------------------------------

"""
Return the factorial ``n! = ∏(1, n)``.
"""
Factorial(n) = Nemo.factorial(ZZ(n))

"""
(SIGNATURES)
"""
Laplace(s, k) = factorial(k) * coeff(s, k)

"""
Return the Laplace expansion of a bivariate exponential
generating function as a power series, or, if 'coeff=true',
as the coefficients of the series.
"""
function EgfExpansion(prec, gf::Function, coeff=true)
    R, x = QPolyRing("x")
    S, t = PowerSeriesRing(R, prec + 1, "t")
    ser = gf(x, t)
    P = Laplace.(ser, 0:prec - 1)
    coeff ? Coefficients.(P) : P
end

"""
Return the Laplace expansion of a bivariate exponential
generating function as a power series.
"""
EgfExpansionCoeff(prec, gf::Function) = EgfExpansion(prec, gf, true)

"""
Return the Laplace expansion of a bivariate exponential generating
function as a list of the coefficients of the power series.
"""
EgfExpansionPoly(prec, gf::Function) = EgfExpansion(prec, gf, false)

function OrthoPoly(dim::Int, s::Function, t::Function)
    T = ZTri(dim, reg=true) # regular triangle
    for n ∈ 1:dim T[n][n] = 1 end
    for n ∈ 2:dim
        u(k) = k == 0 || k == n ? 0 : T[n - 1][k]
        for k ∈ 1:n - 1
            T[n][k] = (u(k - 1) + s(k - 1) * u(k) + t(k) * u(k + 1))
        end
    end
    T
end

"""
(SIGNATURES)
"""
function DelehamΔ(dim::Int, s::Function, t::Function)
    T = ZTri(dim)
    R, x = ZPolyRing("x")
    A = [R(s(k) + x * t(k)) for k ∈ 0:dim - 2]
    C = [R(1) for i ∈ 0:dim]
    C[1] = R(0)
    m = 1
    for k ∈ 0:dim - 1
        for j ∈ k + 1:-1:2
            C[j] = C[j - 1] + C[j + 1] * A[j - 1]
        end
        T[m] = [coeff(C[2], j) for j ∈ 0:k]
        m += 1
    end
    T
end

# --------------------------------------------------

"""
(SIGNATURES)
"""
function RiordanSquare(dim::Int, s::Function)
    T = ZTri(dim, reg=true)
    for n ∈ 0:dim - 1
        T[n + 1][1] = s(n)
    end
    for k ∈ 1:dim - 1, m ∈ k:dim - 1
        T[m + 1][k + 1] = sum(T[j][k] * T[m - j + 2][1] for j ∈ k:m)
    end
    T
end

"""
(SIGNATURES)
"""
function ExpRiordanSquare(dim::Int, s::Function)
    R = RiordanSquare(dim, s)
    u = ZZ(1)
    for k ∈ 1:dim - 1
        u *= k
        for m ∈ 1:k
            j = (m == 1 ? u : div(j, m - 1))
            R[k + 1][m] *= j
        end
    end
    R
end

# --------------------------------------------------

"""
(SIGNATURES)
"""
function check(invM, dim)
    for n = 1:dim
        for k = 1:n
            if denominator(invM[n, k]) != 1
                return false
            end
        end
    end
    return true
end

"""
(SIGNATURES)
"""
function InverseTriangle(T)
    dim = length(T)
    M = zeros(QQ, dim, dim)
    for n in 1:dim, k in 1:n
        M[n, k] = QQ(T[n][k])
    end
    S = MatrixSpace(QQ, dim, dim)

    invM = try
        inv(S(M))
        catch e
        if isa(e, Exception)
            if WARNING_ON_NOTINVERTIBLE
                @warn("\n\n Not invertible!\n\n")
            end
            return ℤSeq[]
        end
    end

    # check if all entries are integer.
    if check(invM, dim) != true
        if WARNING_ON_NOTINVERTIBLE
            @warn("\n\n Inverse not in MatrixSpace(ZZ)!\n\n")
        end
        return ℤSeq[]
    end

    return [[numerator(invM[n, k]) for k = 1:n] for n = 1:dim]
end

"""
Alias for InverseTriangle
"""
Inverse(T::ℤTri) = InverseTriangle(T)

"""
Return the row reversed triangle.
"""
Reverse(T::ℤTri) = reverse.(T)

"""
(SIGNATURES)
"""
function RevInv(T::ℤTri)
    I = Inverse(T)
    return I != [] ? Reverse(I) : ℤSeq[]
end

"""
(SIGNATURES)
"""
InvRev(T::ℤTri) = Inverse(Reverse(T))

"""
(SIGNATURES)
"""
function Flat(T::ℤTri)
    T === [] && return []
    Empty(s) = isempty(s) ? [ZZ(0)] : s
    [z for t ∈ T for z ∈ Empty(t)]
end

"""
(SIGNATURES)
"""
function DeFlat(S::ℤSeq)
    S == [] && return []
    len = length(S)
    n = 1
    while div(n*(n + 1), 2) < len
        n += 1
    end
    T = ZTri(n-1)
    #[z for z ∈ T[k, tria(n)]]
end

"""
Print the array without typeinfo.
"""
function Println(io, v::AbstractVector, newline=true)
    print(io, "[")
    for (i, el) ∈ enumerate(v)
        i > 1 && print(io, ", ")
        print(io, el)
    end
    newline ? println(io, "]") : print(io, "]")
end

# TransformTriangle(Trans::Function, dim::Int, f::Function) =  ℤTri(dim, n -> Trans(n+1, x -> f(x, n)))

# Print the array without typeinfo.
Println(io, str::String) = println(io, str)
Println(T::ℤSeq) = Println(IOContext(stdout), T)
Println(T::ℤTri) = Println.(IOContext(stdout), T)

Print(io, str::String) = print(io, str)
Print(T::ℤSeq) = Println(IOContext(stdout), T, false)
Print(T::ℤTri) = Println.(IOContext(stdout), T, false)

"""
A recursive triangle `RecTriangle` is a subtype of `AbstractTriangle`. The rows
of the triangle are generated by a function `gen(n, k, prevrow)` defined for
``n ≥ 0`` and ``0 ≤ k ≤ n``. The function returns value of type ℤInt.

The parameter prevrow is a function which returns the values of `row(n-1)` of
the triangle and 0 if ``k < 0`` or ``k > n``. The function prevrow is provided
by an instance of `RecTriangle` and must not be defined by the user.
"""
struct RecTriangle <: AbstractTriangle
    dim::Int
    A::ℤSeq
    gen::Function  # generating function
    function RecTriangle(dim::Int, gen::Function)
        dim <= 0 && throw(ArgumentError("dim has to be a positive integer."))
        new(
            dim,
            fill(ZZ(0), dim),
            (n::Int, k::Int, f::Function) -> gen(n, k, f)::ℤInt,
        )
    end
end

function Base.iterate(T::RecTriangle)
    top = T.gen(0, 0, () -> undef)
    T.A[1] = ZZ(top)
    ([top], 1)
end

"""
Return the next row of the triangle.
"""
function Base.iterate(T::RecTriangle, n)
    @inline prevrow(A, n) = (k) -> (k < 0 || k > n) ? ZZ(0) : A[k + 1]
    n == T.dim && return nothing
    F(n, k) = T.gen(n, k, prevrow(T.A, n))
    row = [F(n, k) for k ∈ 0:n]
    for k ∈ 1:n + 1
        T.A[k] = row[k]
    end
    (row, n + 1)
end

Base.length(R::RecTriangle) = R.dim
Base.eltype(R::RecTriangle) = ℤInt

function PolyArray(T::ℤTri)
    P = Polynomial(T)
    # P = ReversePolynomial(T)
    dim = length(T)
    U = ZTri(dim)
    for n = 1:dim
        p = P[n]
        eva = [Evaluate(p, k) for k in 0:dim - 1]
        # Make sure that all denominators are 1.
        # @assert all(d -> d == ZZ(1), denominator.(eva)) "is not an integer"
        U[n] = numerator.(eva)
    end
    U
end

"""
(SIGNATURES)
"""
function PolyTriangle(T::ℤTri)
    A = PolyArray(T)
    U = ZTri(length(T))
    for n = 1:length(T)
        U[n] = [A[n - k][k + 1] for k = 0:n - 1]
    end
    U
end

PolyTri(T::ℤTri) = PolyTriangle(T)

"""
(SIGNATURES)
"""
function PolyValue(T::ℤTri, k::Int)
    P = PolyArray(T)
    [P[n][k + 1] for n = 1:length(P)]
end

"""
(SIGNATURES)
"""
PolyVal2(T::ℤTri) = PolyValue(T, 2)

"""
(SIGNATURES)
"""
PolyVal3(T::ℤTri) = PolyValue(T, 3)

"""
(SIGNATURES)
"""
function PolyValn(T::ℤTri)
    P = PolyArray(T)
    [P[n][n] for n = 1:length(P)]
end

"""
(SIGNATURES)
"""
const TRAITS = Function[
    Flat,
    Reverse,
    Inverse,
    RevInv,
    InvRev,
    DiagTri,
    PolyTri,
    sum,
    EvenSum,
    OddSum,
    AltSum,
    DiagSum,
    Middle,
    Central,
    LeftSide,
    RightSide,
    PosHalf,
    NegHalf,
    PolyVal2,
    PolyVal3,
    PolyValn,
    BinConv,
    IBinConv,
    TransSqrs,
    TransNat0,
    TransNat1
    ]

# START-TEST-########################################################
# using Test

function test()
end

function demo()
    Apery = ℤInt[
        1,
        5,
        73,
        1445,
        33001,
        819005,
        21460825,
        584307365,
    ]

    BinomialTransform(Apery) |> println
    Binomial(Apery) |> println

    typeof(Apery) |> println
    println(typeof(Apery) === ℤSeq)
    typeof(Binomial(Apery)) |> println

    T = ℤSeq[
        ℤInt[1],
        ℤInt[1, 1],
        ℤInt[2, 4, 1],
        ℤInt[6, 18, 9, 1],
        ℤInt[24, 96, 72, 16, 1],
        ℤInt[120, 600, 600, 200, 25, 1],
        ℤInt[720, 4320, 5400, 2400, 450, 36, 1]
    ]

    println("Trans")
    Trans(T, [ZZ(1) for n = 0:33]) |> println

    DiagonalTriangle(T) |> println
    println(ZTri(3, reg=true))
    println(ZTri(3))

    q = QQ(3, 4)
    q |> println
    typeof(q) |> println

    K(n) = [ZZ(1) for _ in 0:n]
    [LinMap(K, K(n)) for n in 0:9] |> println
    L(n) = [ZZ(k) for k in 0:n]
    [LinMap(L, L(n)) for n in 0:9] |> println

    t = Telescope(ℤInt[1,2,3,4,5])
    typeof(t) |> println
    collect(t) |> println

    egfBernoulli(x, t) = divexact(t * exp(x * t), 1 - exp(-t))
    BernoulliPolynomials(n) = EgfExpansionPoly(n, egfBernoulli)
    BernoulliPolynomials(6) |> println

    CatalanTriangle(dim) = OrthoPoly(dim, n -> 0, n -> 1)
    CatalanTriangle(7) |> println
    DiagonalTriangle(CatalanTriangle(9)) |> println

    P = Polynomial(T)
    p = P[3]
    p |> println
    typeof(p) |> println
    s = subst(p, 2)
    s |> println
    s = subst(p, QQ(1, 3))
    s |> println

    T = Telescope(6, n -> ZZ(n))
    println(isa(collect(T), ℤTri))
end

function perf()
end

function main()
    test()
    demo()
    perf()
end

main()

end # module


text:
#=Generator function for TSP that takes the weight Matrix
and returns a function that evaluates the fitness of a single path=#

function fullFit(costs::AbstractMatrix{Float64})
  N = size(costs,1)
  function fFit(tour::Vector{T}) where T<:Integer
    @assert length(tour) == N "Tour must be of length $N"
    @assert isperm(tour) "Not a valid tour, not a permutation"
    #distance = weights[from,to] (from,to) in tour
    distance = costs[tour[N],tour[1]]
    for i in 1:N-1
      @inbounds distance += costs[tour[i],tour[i+1]]
    end
    return distance
  end
  return fFit
end

function partFit(costs::AbstractMatrix{Float64})
  N = size(costs,1)
  function pFit(tour::Vector{T}) where T<:Integer
    n = length(tour)
    #distance = weights[from,to] (from,to) in tour
    distance = n == N ? costs[tour[N],tour[1]] : zero(Float64)
    for i in 1:n-1
      @inbounds distance += costs[tour[i],tour[i+1]]
    end
    return distance
  end
  return pFit
end


text:
using BenchmarkTools, Plots, SNOPT7
# Car escape
T = Float64

# options
max_con_viol = 1.0e-8
verbose=false

opts_ilqr = iLQRSolverOptions{T}(verbose=verbose,
    live_plotting=:off)

opts_al = AugmentedLagrangianSolverOptions{T}(verbose=verbose,
    opts_uncon=opts_ilqr,
    cost_tolerance=1.0e-4,
    cost_tolerance_intermediate=1.0e-2,
    constraint_tolerance=max_con_viol,
    penalty_scaling=50.,
    penalty_initial=10.)

opts_pn = ProjectedNewtonSolverOptions{T}(verbose=verbose,
    feasibility_tolerance=max_con_viol)

opts_altro = ALTROSolverOptions{T}(verbose=verbose,
    opts_al=opts_al,
    opts_pn=opts_pn,
    projected_newton=true,
    projected_newton_tolerance=1.0e-3)

opts_ipopt = DIRCOLSolverOptions{T}(verbose=verbose,
    nlp=Ipopt.Optimizer(),
    feasibility_tolerance=max_con_viol)

opts_snopt = DIRCOLSolverOptions{T}(verbose=verbose,
    nlp=SNOPT7.Optimizer(),
    feasibility_tolerance=max_con_viol)

x0 = Problems.car_3obs.x0
xf = Problems.car_3obs.xf

# ALTRO w/ Newton
prob_altro = copy(Problems.car_3obs)
@time p1, s1 = solve(prob_altro, opts_altro)
@benchmark p1, s1 = solve($prob_altro, $opts_altro)
max_violation_direct(p1)
Problems.plot_car_3obj(p1.X,x0,xf)

# DIRCOL w/ Ipopt
prob_ipopt = copy(Problems.car_3obs)
rollout!(prob_ipopt)
prob_ipopt = update_problem(prob_ipopt,model=Dynamics.car) # get continuous time model)
@time p2, s2 = solve(prob_ipopt, opts_ipopt)
@benchmark p2, s2 = solve($prob_ipopt, $opts_ipopt)
max_violation_direct(p2)
Problems.plot_car_3obj(p2.X,x0,xf)

# DIRCOL w/ SNOPT
prob_snopt = copy(Problems.car_3obs)
rollout!(prob_snopt)
prob_snopt = update_problem(prob_snopt,model=Dynamics.car) # get continuous time model
@time p3, s3 = solve(prob_snopt, opts_snopt)
@benchmark p3, s3 = solve($prob_snopt, $opts_snopt)
max_violation_direct(p3)
Problems.plot_car_3obj(p3.X,x0,xf)
plot(p3.U)

###############################################
#            Create PGF Plot                  #
###############################################

color_obs= "gray"
style = "color=$color_obs, fill=$color_obs"
p = [PGF.Plots.Circle(circle..., style=style) for circle in Problems.circles_3obs]
t1 = trajectory_plot(p1, mark="*", legendentry="ALTRO", style="very thick, color=$col_altro, mark options={fill=$col_altro}");
t2 = trajectory_plot(p2, mark="*", legendentry="Ipopt", style="very thick, color=$col_ipopt, mark options={fill=$col_ipopt}");
t3 = trajectory_plot(p3, mark="*", legendentry="SNOPT", style="very thick, color=$col_snopt, mark options={fill=$col_snopt}");

goal = ([x0[1], xf[1]],
        [x0[2], xf[2]])
z = ["start","end"]
g = PGF.Plots.Scatter(goal[1], goal[2], z,
    scatterClasses="{start={yellow, mark=*, yellow, scale=2},
        end={mark=square*, red, scale=2}}",);

a = Axis([p; t3; t2; t1; g],
    xmin=-0.1, ymin=-1, xmax=1.5, ymax=1,
    axisEqualImage=true,
    legendPos="north west",
    hideAxis=true)

# Save to tikz format
# NOTE: To fix the problem with the legend for the start and goal points, replace \addplot+ with \addplot in the tikz file
save(joinpath(paper,"3obs_traj.tikz"), a, include_preamble=false)

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
using OrderedCollections

abstract type PlottingRecipe end

Plots.plot(pr::PlottingRecipe, args...) = plot!(plot(), pr, args...)

# =============================================================================
# hospitalizations

struct ObservationsPlottingRecipe{Tsd, Ted, To, Te, Tl} <: PlottingRecipe
    startdate ::Tsd
    enddate   ::Ted
    observed  ::To
    expected  ::Te
    label     ::Tl
end

function expected(data, gp, label)
    @assert label in ["cases", "hospitalizations", "deaths"]
    dates  = data.dates
    if label == "cases"
        values = gp.expected_daily_cases
        return (; dates, values )
    elseif label == "hospitalizations"
        values = gp.expected_daily_hospit
        return (; dates, values )
    else
        values = gp.expected_daily_deaths
        return (; dates, values )
    end
end

function observed(data, label)
    @assert label in ["cases", "hospitalizations", "deaths"]
    dates  = data.hospit.date
    if label == "cases"
        values = data.cases.country
        return (; dates, values )
    elseif label == "hospitalizations"
        values = data.hospit.country
        return (; dates, values )
    else
        values = data.deaths.country
        return (; dates, values )
    end
end

function startdate(data::National.Data, label::String)
    @assert label in ["cases", "hospitalizations", "deaths"]
    s = if label == "cases"
        data.turing_data.casemodel.start
    elseif label == "hospitalizations"
        data.turing_data.hospitmodel.start
    else
        data.turing_data.deathmodel.start
    end
    return data.dates[s]
end

function enddate(data::National.Data, label::String)
    @assert label in ["cases", "hospitalizations", "deaths"]
    s = if label == "cases"
        data.turing_data.casemodel.stop
    elseif label == "hospitalizations"
        data.turing_data.hospitmodel.stop
    else
        data.turing_data.deathmodel.stop
    end
    return data.dates[s]
end

function ObservationsPlottingRecipe(data::National.Data, gp, label)
    e = expected(data, gp, label)
    o = observed(data, label)
    sd = startdate(data, label)
    ed = enddate(data, label)
    return ObservationsPlottingRecipe( sd, ed, o, e, label)
end

function Plots.plot!(p::Plots.Plot, r::ObservationsPlottingRecipe)
    o  = r.observed
    e  = r.expected
    ed = r.enddate
    sd = r.startdate

    plot!(p, o.dates, o.values, α=0.5, lc=:match, lab="observed $(r.label)", c=:midnightblue, lw=4, ylab="cases")
    vline!(p, [ed], lab="end observations", lw=2, lc=:black, hover="$ed")
    vline!(p, [sd], lab="start observations", lw=2, lc=:black, hover="$sd", ls=:dash)
    plot_confidence_timeseries!(p, e.dates, e.values; label = "expected $(r.label)") #Dict(hover=>strdates)
end

# ============================================================================
# seropos

struct SeroPlottingRecipe{E, O, L} <: National.PlottingRecipe
    expected  ::E
    observed  ::O
    label     ::L
end

function SeroPlottingRecipe(data::National.Data, gp, label)
    @unpack dates, mean, CIs, delay, population = data.turing_data.seromodel
    observed   = let
        mean *= 100
        error = [(m-l*100, u*100-m) for (m,(l,u)) in zip(mean,CIs)]
        (; dates, mean, error )
    end

    expected = let
        values = gp.expected_seropos  ./ population .* 100
        dates  = data.dates + Day(delay)
        (; dates, values)
    end
    SeroPlottingRecipe( expected, observed, label)
end

function Plots.plot!(p::Plots.Plot, r::SeroPlottingRecipe)
    e = r.expected
    o = r.observed

    plot_confidence_timeseries!(p, e.dates, e.values; r.label)
    y  = o.mean
    ye = o.error
    scatter!(p, o.dates, y, yerror=ye, hover=["$d: $val" for (d, val) in zip(o.dates, y)])
end
# ============================================================================
# reproduction number

struct RtPlottingRecipe{Tlo, Ted, Te, Tl} <: National.PlottingRecipe
    lockdown  ::Tlo
    enddates  ::Ted
    expected  ::Te
    label     ::Tl
end

function RtPlottingRecipe(data::National.Data, gp, label)
    dates      = data.dates
    values     = label == "reproduction number" ? gp.Rt : gp.effective_Rt
    expected   = (; dates, values )
    lockdown   = Date(data.lockdown)
    enddates   = Date(data.observations_end)
    RtPlottingRecipe( lockdown, enddates, expected, label)
end

function Plots.plot!(p::Plots.Plot, r::RtPlottingRecipe)
    e  = r.expected
    ed = r.enddates
    lo = r.lockdown

    vline!(p, [ed], lab="end observations", lw=2, lc=:black, hover="$ed")
    vline!(p, [lo], lab="lockdown", lw=2, lc=:black, hover="$lo", ls=:dash)
    plot_confidence_timeseries!(p, e.dates, e.values; r.label) #Dict(hover=>strdates)
end
# ============================================================================
# infection ascertainment rate

struct IARPlottingRecipe{Tso, Ted, Te} <: National.PlottingRecipe
    startdate ::Tso
    enddates  ::Ted
    expected  ::Te
end

function IARPlottingRecipe(data::National.Data, generated_posterior, args...)
    dates      = data.dates
    values     = generated_posterior.iar
    expected   = (; dates, values )
    startdate  = dates[ data.turing_data.iar_start_idx ]
    enddate    = Date(data.observations_end)
    IARPlottingRecipe( startdate, enddate, expected)
end

function Plots.plot!(p::Plots.Plot, r::IARPlottingRecipe)
    e  = r.expected
    ed = r.enddates
    sd = r.startdate

    vline!(p, [ed], lab="end observations", lw=2, lc=:black, hover="$ed")
    vline!(p, [sd], lab="start case observations", lw=2, lc=:black, hover="$sd", ls=:dash)
    plot_confidence_timeseries!(p, e.dates, e.values; label = "infections ascertainment rate") #Dict(hover=>strdates)
end
# ============================================================================
# region plot

struct OverviewPlottingRecipe{Tr,Tt} <: National.PlottingRecipe
    recipes::Tr
    titles ::Tt
end

posterior2recipe = OrderedDict(
    :expected_daily_cases  => National.ObservationsPlottingRecipe,
    :expected_daily_hospit => National.ObservationsPlottingRecipe,
    :expected_daily_deaths => National.ObservationsPlottingRecipe,
    :Rt                    => National.RtPlottingRecipe,
    :iar                   => National.IARPlottingRecipe,
    :effective_Rt          => National.RtPlottingRecipe,
    :expected_seropos      => National.SeroPlottingRecipe,
)

posterior2label = OrderedDict(
    :expected_daily_cases  => "cases",
    :expected_daily_hospit => "hospitalizations",
    :expected_daily_deaths => "deaths",
    :Rt                    => "reproduction number",
    :iar                   => "infections ascertainment rate",
    :effective_Rt          => "effective reproduction number",
    :expected_seropos      => "sero positive"
)

function OverviewPlottingRecipe(data::National.Data, generated_posterior)
    ks = keys(generated_posterior)
    recipes = Vector{National.PlottingRecipe}()
    labels  = Vector{String}()
    for (k,recipe) in posterior2recipe
        if k in ks
            label = posterior2label[k]
            r     = recipe(data, generated_posterior, label)
            push!(recipes, r)
            push!(labels, label)
        end
    end
    OverviewPlottingRecipe(recipes, labels)
end

function Plots.plot(r::OverviewPlottingRecipe)
    plots  = Vector{Plots.Plot}()
    nplots = length(r.recipes)
    for (recipe, title) in zip(r.recipes, r.titles)
        p = plot(; xaxis = true, legend = :outertopright, title)
        plot!(p, recipe)
        push!(plots, p)
    end
    plot(plots..., layout=(nplots,1), size=(1500, nplots*250), sharex=true, link=:x)
end

Plots.plot(data::National.Data, generated_posterior) =
    plot( OverviewPlottingRecipe(data, generated_posterior) )
# Plots.plot(data::National.Data, generated_posterior) = plot!(plot(), data, generated_posterior)

## ============================================================================
# out dated
function plot_results(data, expected_daily_cases, expected_daily_deaths, expected_daily_hospit, Rt, cumulative_cases, iar)
    l = @layout [a; b; c; d; e; f]

    xlims   = (data.dates[1], data.dates[end])
    # p1 = plot(; xaxis = false, xlims, legend = :outertopright)
    # plot_death_counts!(p1, data)
    p2 = plot(; xaxis = false, xlims, legend = :outertopright)
    plot_deaths_posterior!(p2, data, expected_daily_deaths)
    p3 = plot(; xaxis = false, xlims, legend = :outertopright)
    plot_cases_posterior!(p3, data, expected_daily_cases)
    p4 = plot(; xaxis = false, xlims, legend = :outertopright)
    plot_hospit_posterior!(p4, data, expected_daily_hospit)
    # plot_cumdeaths_posterior!(p5, data, edeaths)
    p5 = plot(; xaxis = true, xlims, legend = :outertopright)
    plot_cumcases_posterior!(p5, data, cumulative_cases)
    p6 = plot(; xaxis = false, xlims, legend = :outertopright)
    plot_Rt_posterior!(p6, data, Rt)
    p7 = plot(; xaxis = false, xlims, legend = :outertopright)
    plot_iar_posterior!(p7, data, iar)
    tickdates = data.dates[1]:Month(1):data.dates[end]
    p = plot(p3, p4, p2, p5, p6, p7, layout=l, size=(1500, 1000), sharex=true, link=:x, xticks = tickdates, legend=false)

    #dates = data.dates[1:20:end] .|> Symbol .|> String
    #Plots.xticks!(1:20:ntot, dates, xrotation=45)
    p
end

function plot_death_counts!(p, data)
    title!(p, "death counts")
    cases_start = data.num_observations - data.turing_data.cases_start
    tstart = data.dates[data.turing_data.epidemic_start]
    tNobs  = data.dates[data.num_observations]
    tcases = data.dates[cases_start]
    plot!(p, data.dates, data.deaths, lc=:match, c=:midnightblue, label=nothing, lw=4) #, label="Observed daily deaths"
    vline!(p, [tstart], label="epidemic start", linewidth=2, lc=:green, hover="$(tstart)")
    vline!(p, [tNobs], label="end of observations", linewidth=2, lc=:black, hover="$(tNobs)")
    vline!(p, [tcases], label="start with case observations", linewidth=2, lc=:black, ls=:dash, hover="$(tcases)")
    # xlims!(p, data.dates[1], data.dates[end])
end

function plot_deaths_posterior!(p, data, expected_daily_deaths)
    edeaths = hcat(expected_daily_deaths...)
    tNobs  = data.dates[data.num_observations]
    title!(p, "posterior death counts")
    plot!(p, data.dates, data.deaths, alpha=0.5, lc=:match, c=:midnightblue, label=nothing, lw=4) #, label="Recorded daily deaths (observed)"
    vline!(p, [tNobs], label="end of observations", linewidth=2, lc=:black, hover="$(tNobs)")
    plot_confidence_timeseries!(p, data.dates, edeaths; label = "Expected daily deaths") #Dict(hover=>strdates)
end

function plot_hospit_posterior!(p, data, expected_daily_hospit)
    ehospit  = hcat(expected_daily_hospit...)
    tNobs  = data.dates[data.num_observations]
    title!(p, "posterior hospitalization counts")
    plot!(p, data.dates, data.hospit, alpha=0.5, lc=:match, label=nothing, c=:midnightblue, lw=4) #, label="Recorded daily hospitalizations (observed)"
    vline!(p, [tNobs], label="end of observations", linewidth=2, lc=:black, hover="$(tNobs)")
    plot_confidence_timeseries!(p, data.dates, ehospit; label = "Expected daily hospitalizations") #Dict(hover=>strdates)
end

function plot_Rt_posterior!(p, data, Rt)
    eRt    = hcat(Rt...)
    tNobs  = data.dates[data.num_observations]
    title!(p, "R_t")
    plot_confidence_timeseries!(p, data.dates, eRt; no_label = true)
    qs = [quantile(v, [0.025, 0.975]) for v in eachrow(eRt)]
    lq, hq = (eachrow(hcat(qs...))..., )
    ylims!(p, 0, maximum(hq) + 0.1)
    vline!(p, [tNobs], label="end of observations", linewidth=2, lc=:black, hover="$(tNobs)")
end

function plot_iar_posterior!(p, data, iar)
    iar    = hcat(iar...)
    size(iar, 1) == 1 && (iar = repeat(iar, outer=( length(data.dates), 1)) )
    tNobs  = data.dates[data.num_observations]
    title!(p, "infection ascertainment rate")
    plot_confidence_timeseries!(p, data.dates, iar; no_label = true)
    qs = [quantile(v, [0.025, 0.975]) for v in eachrow(iar)]
    lq, hq = (eachrow(hcat(qs...))..., )
    ylims!(p, 0, maximum(hq) + 0.1)
    vline!(p, [tNobs], label="end of observations", linewidth=2, lc=:black, hover="$(tNobs)")
end

function plot_cases_posterior!(p, data, expected_daily_cases)
    ecases  = hcat(expected_daily_cases...)
    cases_start = data.num_observations - data.turing_data.num_case_obs
    tcases = data.dates[cases_start]
    tNobs  = data.dates[data.num_observations]
    title!(p, "posterior case counts")
    plot!(p, data.dates, data.cases, alpha=0.5, lc=:match, label=nothing, c=:midnightblue, lw=4)
    vline!(p, [tcases], label="start with case observations", linewidth=2, lc=:black, ls=:dash, hover="$(tcases)")
    vline!(p, [tNobs], label="end of observations", linewidth=2, lc=:black, hover="$(tNobs)")
    # plot!(p, [d1, d2, d2, d1], [low, low, up, up], fc=:orange, fα=0.4, leg=nothing, hover="$d1 - $d2", seriestype=:shape, lc=:orange)
    plot_confidence_timeseries!(p,data.dates, ecases; label = "observed daily cases")
    # plot_confidence_timeseries!(p,data.dates, ecases; label = "expected daily cases")
end

function plot_cumdeaths_posterior!(p, data, expected_daily_deaths)
    edeaths = hcat(expected_daily_deaths...)
    tNobs  = data.dates[data.num_observations]
    title!(p, "posterior cumulative death counts")
    vline!(p, [tNobs], label="end of observations", linewidth=2, lc=:black, hover="$(tNobs)")
    plot_confidence_timeseries!(p, data.dates, cumsum(edeaths; dims=1); label = "Expected cumulative deaths")
    Plots.plot!(p, data.dates, cumsum(data.deaths; dims=1), label="Recorded cumulative deaths", color=:red)
end

function plot_cumcases_posterior!(p, data, cumulative_cases)
    ecumul  = hcat(cumulative_cases...)
    tNobs  = data.dates[data.num_observations]
    title!(p, "posterior cumulative cases counts")
    vline!(p, [tNobs], label="end of observations", linewidth=2, lc=:black, hover="$(tNobs)")
    plot_confidence_timeseries!(p, data.dates, ecumul ./ data.pop .* 100; label = "Expected cumulative cases")
    # plot_confidence_timeseries!(p, data.dates, cumsum(ecases, dims=1) ./ data.pop .* 100; label = "test")
    Plots.plot!(p, data.dates, cumsum(data.cases; dims=1) ./ data.pop .* 100, label="Recorded cumulative cases", color=:red)
    plot_seroprev!(p, data)

    # plot_seroprev!(p, Date("2020-06-09"), Date("2020-06-10"), 1.2, .7, 1.7)
    # plot_seroprev!(p, Date("2020")+Week(34), Date("2020")+Week(36), 2.2, 1.8, 2.6)
    # plot_seroprev!(p, Date("2020")+Week(38), Date("2020")+Week(51), 3.9, 3.3, 4.6)
    #TODO include seroprevalence study: 34 (1.2%) positive with 95% CI: 0,7%-1,7% in week 20-21
    #TODO include seroprevalence study: 2.2% positive with 95% CI: 1,8% - 2,6% in week 34-36
    #TODO include seroprevalence study: 3.9% positive with 95% CI: 3,3% - 4,6% in week 38-51
end

function plot_seroprev!(p, data)
    seroprev = data.seroprev
    le = [s.country-s.CI[1] for s in eachrow(seroprev)] * 100
    ue = [s.CI[2]-s.country for s in eachrow(seroprev)] * 100
    scatter!(p, seroprev.date, seroprev.country * 100, yerror=(le,ue), hover=["$d" for d in seroprev.date])
end

function plot_seroprev!(p, d1, d2, med, up, low)
    # rec = [(d1, low), (d1, low), (d2, up), (d2, up)]
    plot!(p, [d1, d2, d2, d1], [low, low, up, up], fc=:orange, fα=0.4, leg=nothing, hover="$d1 - $d2", seriestype=:shape, lc=:orange)
    plot!(p, [d1,d2], [med, med], lc=:orange, leg=nothing, lw=4, hover="$med")
end

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
#!/usr/bin/env julia6
using ReadGlobal, MyStats, JLD2

function histstdindices(input::AbstractString,nbins::Integer,nstd::Number)
    field = readfield(input)
    med,stdv,_,_ = read_info(input*".info")
    minf = med - nstd*stdv
    maxf = med + nstd*stdv
    ind = hist_indices(field,minf,maxf,nbins)
    bins = Bins(minf,maxf,nbins)
    @save "pdf_std$(nstd)_indices_$(input)_nbins$(nbins).jld2" ind bins
end

function main(f,ns,nstds)
    n = parse(Int,ns)
    nstd = parse(Float64,nstds)
    histstdindices(f,n,nstd)
end

function main(f,ns)
    n = parse(Int,ns)
    histstdindices(f,n,2.)
end

main(f) = histstdindices(f,50,2.)

main(ARGS...)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct DiskUpdateProperties <: SwaggerModel
    osType::Any # spec type: Union{ Nothing, String } # spec name: osType
    diskSizeGB::Any # spec type: Union{ Nothing, Int32 } # spec name: diskSizeGB
    encryptionSettings::Any # spec type: Union{ Nothing, EncryptionSettings } # spec name: encryptionSettings

    function DiskUpdateProperties(;osType=nothing, diskSizeGB=nothing, encryptionSettings=nothing)
        o = new()
        validate_property(DiskUpdateProperties, Symbol("osType"), osType)
        setfield!(o, Symbol("osType"), osType)
        validate_property(DiskUpdateProperties, Symbol("diskSizeGB"), diskSizeGB)
        setfield!(o, Symbol("diskSizeGB"), diskSizeGB)
        validate_property(DiskUpdateProperties, Symbol("encryptionSettings"), encryptionSettings)
        setfield!(o, Symbol("encryptionSettings"), encryptionSettings)
        o
    end
end # type DiskUpdateProperties

const _property_map_DiskUpdateProperties = Dict{Symbol,Symbol}(Symbol("osType")=>Symbol("osType"), Symbol("diskSizeGB")=>Symbol("diskSizeGB"), Symbol("encryptionSettings")=>Symbol("encryptionSettings"))
const _property_types_DiskUpdateProperties = Dict{Symbol,String}(Symbol("osType")=>"String", Symbol("diskSizeGB")=>"Int32", Symbol("encryptionSettings")=>"EncryptionSettings")
Base.propertynames(::Type{ DiskUpdateProperties }) = collect(keys(_property_map_DiskUpdateProperties))
Swagger.property_type(::Type{ DiskUpdateProperties }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_DiskUpdateProperties[name]))}
Swagger.field_name(::Type{ DiskUpdateProperties }, property_name::Symbol) =  _property_map_DiskUpdateProperties[property_name]

const _allowed_DiskUpdateProperties_osType = ["Windows", "Linux"]

function check_required(o::DiskUpdateProperties)
    true
end

function validate_property(::Type{ DiskUpdateProperties }, name::Symbol, val)
    if name === Symbol("osType")
        Swagger.validate_param(name, "DiskUpdateProperties", :enum, val, _allowed_DiskUpdateProperties_osType)
    end
end


text: # TODO: Will do it in this week