Enable dynamic parameters

Circuit/Analysis.cs

@@ -238,7 +238,11 @@ public Expression AddUnknownEqualTo(string Name, Expression Eq)
         /// </summary>
         /// <param name="Eq"></param>
         /// <returns></returns>
-        public Expression AddUnknownEqualTo(Expression Eq) { return AddUnknownEqualTo(AnonymousName(), Eq); }
+        public Expression AddUnknownEqualTo(Expression x) {
+            if (x is Constant) return x;
+            if (x is Call c && c.Target is UnknownFunction) return x;
+            return AddUnknownEqualTo(AnonymousName(), x); 
+        }
 
         /// <summary>
         /// Add initial conditions to the system.
@@ -253,6 +257,55 @@ public Expression AddUnknownEqualTo(string Name, Expression Eq)
         /// <returns></returns>
         public string AnonymousName() { return context.AnonymousName(); }
 
+        /// <summary>
+        /// Describes a parameter for the circuit.
+        /// </summary>
+        public class Parameter
+        {
+            private Component of;
+            /// <summary>
+            /// Component the parameter affects.
+            /// </summary>
+            public Component Of { get { return of; } }
+
+            private Expression expr;
+            /// <summary>
+            /// Expression for the parameter in the system of equations for this analysis.
+            /// </summary>
+            public Expression Expression { get { return expr; } }
+
+            Func<double> getter;
+            /// <summary>
+            /// Current value of the parameter.
+            /// </summary>
+            public double Value { get { return getter(); } }
+
+            public Parameter(Component Of, Expression Expression, Func<double> Getter)
+            {
+                of = Of;
+                expr = Expression;
+                getter = Getter;
+            }
+        }
+
+        private List<Parameter> parameters = new List<Parameter>();
+        /// <summary>
+        /// Enumerates the parameters in this analysis.
+        /// </summary>
+        public IEnumerable<Parameter> Parameters { get { return parameters; } }
+
+        /// <summary>
+        /// Create an expression for a parameter in the current context.
+        /// </summary>
+        /// <param name="Name"></param>
+        /// <returns></returns>
+        public Expression AddParameter(Component Of, string Name, Func<double> Getter)
+        {
+            Expression expr = Variable.New(context.Prefix + Name);
+            parameters.Add(new Parameter(Of, expr, Getter)); 
+            return expr;
+        }
+
         private void AddKcl(Dictionary<Expression, Expression> kcl, Expression V, Expression i)
         {
             if (kcl.TryGetValue(V, out var sumi))

Circuit/Component.cs

@@ -21,7 +21,7 @@ public interface IGroupableComponent
     /// <summary>
     /// Interface for components that expose a pot controlled value.
     /// </summary>
-    public interface IPotControl: IGroupableComponent
+    public interface IPotControl : IGroupableComponent
     {
         /// <summary>
         /// The value of the pot.
@@ -32,7 +32,7 @@ public interface IPotControl: IGroupableComponent
     /// <summary>
     /// Interface for components that expose a button controlled value.
     /// </summary>
-    public interface IButtonControl: IGroupableComponent
+    public interface IButtonControl : IGroupableComponent
     {
         void Click();
         int Position { get; set; }

Circuit/Components/Capacitor.cs

@@ -25,7 +25,7 @@ public static Expression Analyze(Analysis Mna, string Name, Node Anode, Node Cat
             V = Mna.AddUnknownEqualTo("V" + Name, V);
             // i = C*dV/dt
             Expression i = C * D(V, t);
-            //i = Mna.AddUnknownEqualTo("i" + Name, i);
+            i = Mna.AddUnknownEqualTo("i" + Name, i);
             Mna.AddPassiveComponent(Anode, Cathode, i);
             return i;
         }

Circuit/Components/Potentiometer.cs

@@ -63,7 +63,7 @@ public void ConnectTo(Node A, Node C, Node W)
 
         public override void Analyze(Analysis Mna)
         {
-            Expression P = VariableResistor.AdjustWipe(wipe, sweep);
+            Expression P = Mna.AddParameter(this, Name, () => VariableResistor.AdjustWipe(Wipe, Sweep));
 
             Expression R1 = Resistance * P;
             Expression R2 = Resistance * (1 - P);

Circuit/Components/VacuumTubes/Triode.cs

@@ -146,7 +146,8 @@ public override void Analyze(Analysis Mna)
                     break;
                 case TriodeModel.Koren:
                     Expression E1 = Ln1Exp(Kp * (1.0 / Mu + Vgk * Binary.Power(Kvb + Vpk * Vpk, -0.5))) * Vpk / Kp;
-                    ip = Mna.AddUnknownEqualTo(Call.If(E1 > 0, 2d * (E1 ^ Ex) / Kg, 0));
+                    ip = Call.If(E1 > 0, 2d * (E1 ^ Ex) / Kg, 0);
+                    ip = Mna.AddUnknownEqualTo("Ip", ip);
 
                     var vg = (Real)Vg;
                     var knee = (Real)Kn;
@@ -156,7 +157,8 @@ public override void Analyze(Analysis Mna)
                     var b = (knee - vg) / (2 * knee * rg1);
                     var c = (-a * Binary.Power(vg - knee, 2)) - (b * (vg - knee));
 
-                    ig = Mna.AddUnknownEqualTo(Call.If(Vgk < vg - knee, 0, Call.If(Vgk > vg + knee, (Vgk - vg) / rg1, a * Vgk * Vgk + b * Vgk + c)));
+                    ig = Call.If(Vgk < vg - knee, 0, Call.If(Vgk > vg + knee, (Vgk - vg) / rg1, a * Vgk * Vgk + b * Vgk + c));
+                    ig = Mna.AddUnknownEqualTo("Ig", ig);
                     ik = -(ip + ig);
                     break;
                 case TriodeModel.DempwolfZolzer:

Circuit/Components/VariableResistor.cs

@@ -64,7 +64,7 @@ public class VariableResistor : TwoTerminal, IPotControl
 
         protected double wipe = 0.5;
         [Serialize, Description("Position of the wiper on this variable resistor, between 0 and 1.")]
-        public double Wipe { get { return wipe; } set { wipe = value; NotifyChanged(nameof(Wipe)); } }
+        public double Wipe { get { return wipe; } set { wipe = value; NotifyChanged(nameof(Wipe)); NotifyChanged(nameof(IPotControl.PotValue)); } }
         // IPotControl
         double IPotControl.PotValue { get { return Wipe; } set { Wipe = value; } }
 
@@ -80,7 +80,7 @@ public class VariableResistor : TwoTerminal, IPotControl
 
         public override void Analyze(Analysis Mna)
         {
-            Expression P = AdjustWipe(wipe, sweep);
+            Expression P = Mna.AddParameter(this, Name, () => AdjustWipe(Wipe, Sweep));
 
             Resistor.Analyze(Mna, Name, Anode, Cathode, (Expression)Resistance * P);
         }

Circuit/Components/VoltageSource.cs

@@ -26,6 +26,9 @@ public static void Analyze(Analysis Mna, string Name, Node Anode, Node Cathode,
         {
             // Unknown current.
             Mna.AddPassiveComponent(Anode, Cathode, Mna.AddUnknown("i" + Name));
+
+            V = Mna.AddUnknownEqualTo(V);
+
             // Set the voltage.
             Mna.AddEquation(Anode.V - Cathode.V, V);
         }
@@ -35,7 +38,10 @@ public static void Analyze(Analysis Mna, string Name, Node Anode, Node Cathode,
             Mna.AddInitialConditions(InitialConditions);
         }
         public static void Analyze(Analysis Mna, Node Anode, Node Cathode, Expression V) { Analyze(Mna, Mna.AnonymousName(), Anode, Cathode, V); }
-        public static void Analyze(Analysis Mna, Node Anode, Node Cathode, Expression V, Arrow InitialConditions) { Analyze(Mna, Mna.AnonymousName(), Anode, Cathode, V, InitialConditions); }
+        public static void Analyze(Analysis Mna, Node Anode, Node Cathode, Expression V, Arrow InitialConditions) 
+        { 
+            Analyze(Mna, Mna.AnonymousName(), Anode, Cathode, V, InitialConditions); 
+        }
 
         public override void Analyze(Analysis Mna) { Analyze(Mna, Name, Anode, Cathode, Voltage); }
 

Circuit/Simulation/Simulation.cs

@@ -6,6 +6,7 @@
 using System.Numerics;
 using System.Reflection;
 using Util;
+using static Circuit.Analysis;
 using LinqExpr = System.Linq.Expressions.Expression;
 using ParamExpr = System.Linq.Expressions.ParameterExpression;
 
@@ -98,6 +99,12 @@ public TransientSolution Solution
         /// </summary>
         public IEnumerable<Expression> Output { get { return output; } set { output = value.ToArray(); InvalidateProcess(); } }
 
+        private Expression[] arguments = new Expression[] { };
+        /// <summary>
+        /// Expressions for parameters.
+        /// </summary>
+        public IEnumerable<Expression> Arguments { get { return arguments; } set { arguments = value.ToArray(); InvalidateProcess(); } }
+
         // Stores any global state in the simulation (previous state values, mostly).
         private Dictionary<Expression, GlobalExpr<double>> globals = new Dictionary<Expression, GlobalExpr<double>>();
         // Add a new global and set it to 0 if it didn't already exist.
@@ -153,7 +160,8 @@ public Simulation(TransientSolution Solution)
         /// <param name="N">Number of samples to process.</param>
         /// <param name="Input">Buffers that describe the input samples.</param>
         /// <param name="Output">Buffers to receive output samples.</param>
-        public void Run(int N, IEnumerable<double[]> Input, IEnumerable<double[]> Output)
+        /// <param name="Parameters">Scalar values for the arguments.</param>
+        public void Run(int N, IEnumerable<double[]> Input, IEnumerable<double[]> Output, IEnumerable<double> Parameters)
         {
             if (_process == null)
                 _process = DefineProcess();
@@ -162,7 +170,7 @@ public void Run(int N, IEnumerable<double[]> Input, IEnumerable<double[]> Output
             {
                 try
                 {
-                    _process(N, n*TimeStep, Input.AsArray(), Output.AsArray());
+                    _process(N, n*TimeStep, Input.AsArray(), Output.AsArray(), Parameters.AsArray());
                     n += N;
                 }
                 catch (TargetInvocationException Ex)
@@ -175,11 +183,11 @@ public void Run(int N, IEnumerable<double[]> Input, IEnumerable<double[]> Output
                 throw new SimulationDiverged("Simulation diverged near t = " + Quantity.ToString(Time, Units.s) + " + " + Ex.At, n + Ex.At);
             }
         }
-        public void Run(int N, IEnumerable<double[]> Output) { Run(N, new double[][] { }, Output); }
-        public void Run(double[] Input, IEnumerable<double[]> Output) { Run(Input.Length, new[] { Input }, Output); }
-        public void Run(double[] Input, double[] Output) { Run(Input.Length, new[] { Input }, new[] { Output }); }
+        public void Run(int N, IEnumerable<double[]> Output, IEnumerable<double> Parameters) { Run(N, new double[][] { }, Output, Parameters); }
+        public void Run(double[] Input, IEnumerable<double[]> Output, IEnumerable<double> Parameters) { Run(Input.Length, new[] { Input }, Output, Parameters); }
+        public void Run(double[] Input, double[] Output, IEnumerable<double> Parameters) { Run(Input.Length, new[] { Input }, new[] { Output }, Parameters); }
 
-        private Action<int, double, double[][], double[][]> _process;
+        private Action<int, double, double[][], double[][], double[]> _process;
         // Force rebuilding of the process function.
         private void InvalidateProcess()
         {
@@ -189,11 +197,12 @@ private void InvalidateProcess()
         // The resulting lambda processes N samples, using buffers provided for Input and Output:
         //  void Process(int N, double t0, double T, double[] Input0 ..., double[] Output0 ...)
         //  { ... }
-        private Action<int, double, double[][], double[][]> DefineProcess()
+        private Action<int, double, double[][], double[][], double[]> DefineProcess()
         {
             // Map expressions to identifiers in the syntax tree.
             var inputs = new List<KeyValuePair<Expression, LinqExpr>>();
             var outputs = new List<KeyValuePair<Expression, LinqExpr>>();
+            var parameters = new List<KeyValuePair<Expression, LinqExpr>>();
 
             // Lambda code generator.
             CodeGen code = new CodeGen();
@@ -203,6 +212,7 @@ private Action<int, double, double[][], double[][]> DefineProcess()
             ParamExpr t = code.Decl(Scope.Parameter, Simulation.t);
             var ins = code.Decl<double[][]>(Scope.Parameter, "ins");
             var outs = code.Decl<double[][]>(Scope.Parameter, "outs");
+            var parms = code.Decl<double[]>(Scope.Parameter, "params");
 
             // Create buffer parameters for each input...
             for (int i = 0; i < input.Length; i++)
@@ -216,6 +226,11 @@ private Action<int, double, double[][], double[][]> DefineProcess()
                 outputs.Add(new KeyValuePair<Expression, LinqExpr>(output[i], LinqExpr.ArrayAccess(outs, LinqExpr.Constant(i))));
             }
 
+            for (int i = 0; i < arguments.Length; i++)
+            {
+                parameters.Add(new KeyValuePair<Expression, LinqExpr>(arguments[i], LinqExpr.ArrayAccess(parms, LinqExpr.Constant(i))));
+            }
+
             Arrow t_t1 = Arrow.New(Simulation.t, Simulation.t - Solution.TimeStep);
 
             // Create globals to store previous values of inputs.
@@ -240,6 +255,9 @@ private Action<int, double, double[][], double[][]> DefineProcess()
             foreach (KeyValuePair<Expression, LinqExpr> i in inputs)
                 code.DeclInit(i.Key, code[i.Key.Evaluate(t_t1)]);
 
+            foreach (KeyValuePair<Expression, LinqExpr> i in parameters)
+                code.DeclInit(i.Key, i.Value);
+
             // Create arrays for linear systems.
             int M = Solution.Solutions.OfType<NewtonIteration>().Max(i => i.Equations.Count(), 0);
             int N = Solution.Solutions.OfType<NewtonIteration>().Max(i => i.UnknownDeltas.Count(), 0);
@@ -409,7 +427,7 @@ private Action<int, double, double[][], double[][]> DefineProcess()
             foreach (KeyValuePair<Expression, GlobalExpr<double>> i in globals)
                 code.Add(LinqExpr.Assign(i.Value, code[i.Key]));
 
-            var lambda = code.Build<Action<int, double, double[][], double[][]>>();
+            var lambda = code.Build<Action<int, double, double[][], double[][], double[]>>();
             return lambda.Compile();
         }
 

Circuit/Simulation/TransientSolution.cs

@@ -79,10 +79,10 @@ public static TransientSolution Solve(Analysis Analysis, Expression TimeStep, IE
             // Define T = step size.
             DynamicNamespace globals = new DynamicNamespace();
             globals.Add("T", h);
-            // Define d[t] = delta function.
+            // Define dq = delta function.
             // TODO: This should probably be centered around 0, and also have an integral of 1 (i.e. a height of 1 / h).
             globals.Add(ExprFunction.New("d", Call.If((0 <= t) & (t < h), 1, 0), t));
-            // Define u[t] = step function.
+            // Define uq = step function.
             globals.Add(ExprFunction.New("u", Call.If(t >= 0, 1, 0), t));
             mna = mna.Resolve(Analysis).Resolve(globals).OfType<Equal>().ToList();
 
@@ -101,6 +101,8 @@ public static TransientSolution Solve(Analysis Analysis, Expression TimeStep, IE
                 .Substitute(dy_dt.Select(i => Arrow.New(i, 0)).Append(Arrow.New(t, 0), Arrow.New(T, 0), SinglePoleSwitch.IncludeOpen))
                 // Use the initial conditions from analysis.
                 .Substitute(Analysis.InitialConditions)
+                // Use the default parameter values.
+                .Substitute(Analysis.Parameters.Select(i => Arrow.New(i.Expression, i.Value)))
                 // Evaluate variables at t=0.
                 .OfType<Equal>(), y.Select(j => j.Substitute(t, 0)));
 
@@ -149,12 +151,13 @@ public static TransientSolution Solve(Analysis Analysis, Expression TimeStep, IE
                 Log.WriteLine(MessageType.Verbose, "Partition unknowns: {0}", String.Join(", ", F.Unknowns));
                 // Find linear solutions for y. Linear systems should be completely solved here.
                 F.RowReduce();
-                IEnumerable<Arrow> linear = F.Solve();
-                if (linear.Any())
+                List<Arrow> linearFwd = new List<Arrow>();
+                List<Arrow> linearBack = new List<Arrow>();
+                F.PartialSolve(linearFwd, linearBack);
+                if (linearFwd.Any())
                 {
-                    linear = Factor(linear);
-                    solutions.Add(new LinearSolutions(linear));
-                    LogExpressions(Log, MessageType.Verbose, "Linear solutions:", linear);
+                    solutions.Add(new LinearSolutions(linearFwd));
+                    LogExpressions(Log, MessageType.Verbose, "Linear solutions:", linearFwd);
                 }
 
                 // If there are any variables left, there are some non-linear equations requiring numerical techniques to solve.
@@ -175,20 +178,23 @@ public static TransientSolution Solve(Analysis Analysis, Expression TimeStep, IE
                     // Find linear solutions for dy. 
                     nonlinear.RowReduce(ly);
                     IEnumerable<Arrow> solved = nonlinear.Solve(ly);
-                    solved = Factor(solved);
 
                     // Initial guess for y[t] = y[t - h].
                     IEnumerable<Arrow> guess = F.Unknowns.Select(i => Arrow.New(i, i.Substitute(t, t - h))).ToList();
-                    guess = Factor(guess);
 
                     // Newton system equations.
                     IEnumerable<LinearCombination> equations = nonlinear.Equations.Buffer();
-                    equations = Factor(equations);
 
                     solutions.Add(new NewtonIteration(solved, equations, nonlinear.Unknowns, guess));
                     LogExpressions(Log, MessageType.Verbose, String.Format("Non-linear Newton's method updates ({0}):", String.Join(", ", nonlinear.Unknowns)), equations.Select(i => Equal.New(i, 0)));
                     LogExpressions(Log, MessageType.Verbose, "Linear Newton's method updates:", solved);
                 }
+
+                if (linearBack.Any())
+                {
+                    solutions.Add(new LinearSolutions(linearBack));
+                    LogExpressions(Log, MessageType.Verbose, "Linear solutions:", linearBack);
+                }
             }
 
             Log.WriteLine(MessageType.Info, "System solved, {0} solution sets for {1} unknowns.",
@@ -206,9 +212,6 @@ public static TransientSolution Solve(Analysis Analysis, Expression TimeStep, IE
         public static TransientSolution Solve(Analysis Analysis, Expression TimeStep, ILog Log) { return Solve(Analysis, TimeStep, new Arrow[] { }, Log); }
         public static TransientSolution Solve(Analysis Analysis, Expression TimeStep) { return Solve(Analysis, TimeStep, new Arrow[] { }, new NullLog()); }
 
-        private static IEnumerable<Arrow> Factor(IEnumerable<Arrow> x) { return x.Select(i => Arrow.New(i.Left, i.Right.Factor())).Buffer(); }
-        private static IEnumerable<LinearCombination> Factor(IEnumerable<LinearCombination> x) { return x.Select(i => LinearCombination.New(i.Select(j => new KeyValuePair<Expression, Expression>(j.Key, j.Value.Factor())))).Buffer(); }
-
         // Shorthand for df/dx.
         protected static Expression D(Expression f, Expression x) { return Call.D(f, x); }