Calculating the std deviation and mean along the y axis on a family of waveform

I can't quite see what the sweep variable is in the legend, but if you plot swapSweep(inlEndpoint(Iout_sink) "mcp") does it plot the 10 points across the x-axis and separate traces for each vcode? If so, you should be able to then do stddev() or average() over the result to get what you want.

I'm assuming inlEndpoint is some custom function anyway - I don't think it's something we provide as far as I can tell.

Andrew.

Hi Andrew:

Sorry, maybe I should give more background. I'm new to skill programming so please bear with me if I make any mistake.

I wrote the inlEndpoint with the help on examples in the custom skill library. The code is shown here.

(defglobalfun inlEndpoint (wave)
        (cond
             ((drIsWaveform wave)
             /* get the x,y axis vector of the waveform*/
             (setq xVec (drGetWaveformXVec wave))
             (setq yVec (drGetWaveformYVec wave))
             /*get length of the vector, error shows up if I put extra parenthesis at the outside?*/
             setq(len (drVectorLength xVec))

             /*get first element*/
             setq(firstEle (drGetElem yVec 0))
             setq(adjustedOffset firstEle)
             /*get last element*/
             setq(lastEle (drGetElem yVec len-1))
             adjustedLsb = (lastEle - firstEle) / (len - 1)

             (setq newY (drCreateVec (drGetWaveformYType wave) len))

             (for pos 0 (sub1 len)
                 (drAddElem newY ( drGetElem(yVec pos) - (adjustedOffset + adjustedLsb * pos)) / adjustedLsb )
             )
             (drCreateWaveform xVec newY)
         )
        ((famIsFamily wave)
             (famMap 'inlEndpoint wave)
             )
             (t
             (error "inlEndpoint: cannot handle %L\n" wave)
             )
)
)

I'm running monte carlo run with multiple corners and the setup is shown in the figure

When I plot the inlEndpoint(iout_sink) it's like the figure in the main post. 

I tried the swapSweep(inlEndpoint(Iout_sink) "mcparamset" 10), and I noticed it changes the x axis to mcparamset, but I want to maintain the x axis the same.

What i'm trying to do is I have 4 corner runs with 10 monte carlo samples, where I have a total of 40 waveforms (10 waveform for a single corner run)

In the end I want to combine the 10 waveform for a single corner into 1 waveform, so total I want to output 4 waveforms.

I'm trying to use a single corner for example. I want to do calculations along the y axis for the mean and variance.

waveform0Yvec [0a,0b,0c]

waveform1Yvec [1a,1b,1c]

waveform2Yvec [2a,2b,2c]

..

...

waveform9Yvec [9a,9b,9c]

-> output to 

waveformOutputYVec[mean of (0a, 1a...9a), mean of (0b,1b...9b ), mean of (0c, 1c...2c)] while mantaining the x-axis

I hope i'm clear in my explanations what i'm trying to do... I don't know if there is any other ways to do this 

The pesudo code I can think of is the following (Not sure if it helps...)

Function definition : OutputMeanAndVaiance (takes familywaveform with multiple corners(4) and multiple monte carlo runs(10)) -> output a familywaveform with multiple corners(4)

for every corner {

    take the corresponding 10 monte carlo waveform

             create meanSigmaYVec with length of the XVec of waveform

             for(index = 0; i < xVector.length; i++) {

                    sum = 0;

                    square = 0;

                    for(monteNumber = 0; monteNumber < 10; pos++) {

                             sum += waveformof[monteNumber]Yvec[index];

                             square += waveformof[monteNumber]Yvec[index] ^ 2;

                    }//end monteNumber

              mean = sum / 10;

              sigma = square/10 - mean^2;

              meanSigmaYVec[index] = mean + sigma;

              }//end traversing the xVec

             plot a waveform with meanSigmaYVec and original XVec

}//end for every corner

Thanks

I know that you don't want to have mcparamset as the x-axis in your plots, but the reason for my suggestion is that if you have a set of curves, then the average or stddev functions will give you the mean or standard deviation across the x-axis. You'll then get one value for each y curve, so I think you'll then get a curve of the mean across the mcparamset with the x-axis being vcode again. Hopefully if all is OK, it will also give you one of these curves for each corner. I don't have something similar to you to test with, but I think it should work.

So I don't believe any custom code is needed. The swapSweep is just to allow you to change around the axes so that you can pick which is the x-axis to find the mean or standard deviation across, and then that number will be plotted against the next outer sweep variable. Worst case you might need to swapSweep the result of the average/stddev to change the axis back to vcode again (if it happens to plot versus corner instead, say), but hopefully that's not necessary.

If this doesn't make sense or you can't make it work, please contact customer support. Point the application engineer to this thread if you do that, and mention that I'd be happy to help on this (by the way, many fo the examples in the calculator function code library come from me!)

Regards,

Andrew

The SwapSweep function from the calculator works fine to produce a waveform from values of all sweep elements, however it only processes a single x-value; it does not fully transpose sweep and x-values.

As I assume that you want the statistics at all available x-values, you can try the attached script. It should be able to produce the waveforms you want. It is however only a fast hack, not really elegant and might also not work on every input. If you do not need the third and forth moments, you can also remove alot.
Call from calculator with:
icFamStat(<option> <list of waveforms or families>)

You can fill in arbitrary waveforms and families, the first processed wave will determine the x-values and unit of the y-values.

The function can be registered with the calculator using the fx button in the function panel (IC617 and later). It can then be found under "Custom function". A GUI builder is at the end, so the function will be available automtically when loaded appropriately.

Example output:

;====================================================================
;; icFamStatCore
;; collect non-central moments of 1st..4th order
;====================================================================
procedure( icFamStatCore(wave work)
	let( (xVec yVec len y1 y2 y3 y4 N v p s)

		xVec = drGetWaveformXVec(wave)
		yVec = drGetWaveformYVec(wave)
		len  = drVectorLength(xVec)

		if( !work->xVec then
			;; initialize
						
			y1= drCreateVec('double len)
			y2= drCreateVec('double len)
			y3= drCreateVec('double len)
			y4= drCreateVec('double len)
						
			for(i 0 len-1
				v = drGetElem(yVec i)

				p = v
				drAddElem(y1 p)
				p = p*v
				drAddElem(y2 p)
				p = p*v
				drAddElem(y3 p)
				p = p*v
				drAddElem(y4 p)
			)
			y1->units  = yVec->units
			y1->name   = "mean"
			y2->units  = yVec->units
			y2->name   = "sqr"
			y3->units  = yVec->units
			y3->name   = "cub"
			y4->units  = yVec->units
			y4->name   = "quad"

			work->xVec = xVec
			work->y1   = y1
			work->y2   = y2
			work->y3   = y3
			work->y4   = y4
			work->N    = makeVector(len 1)
			work->len  = len
		else
			when( xVec!=work->xVec
				error("icFamStat: different x-axis in %L\n" wave)
			)

			for(i 0 len-1
				v = drGetElem(yVec i)
				when(v 
					N = work->N[i] + 1
					work->N[i] = N

					p = v
					;; recurrence formulae (mean of y)
					s = drGetElem(work->y1 i)
					s = s + (p - s) / float(N)
					drSetElem(work->y1 i s)

					p = p*v
					;; recurrence formulae (mean of y^2)
					s = drGetElem(work->y2 i)
					s = s + (p - s) / float(N)
					drSetElem(work->y2 i s)

					p = p*v
					;; recurrence formulae (mean of y^3)
					s = drGetElem(work->y3 i)
					s = s + (p - s) / float(N)
					drSetElem(work->y3 i s)

					p = p*v
					;; recurrence formulae (mean of y^4)
					s = drGetElem(work->y4 i)
					s = s + (p - s) / float(N)
					drSetElem(work->y4 i s)
				)
			)
		); if
	); let
)
;====================================================================

;====================================================================
;; icFamStatHop
;; iterate over families, call core with single waveforms
;; collect data in work
;====================================================================
procedure( icFamStatHop(wave work)
	let( (wf)
		cond(
			( drIsWaveform(wave)
				 icFamStatCore(wave work)
			)
			( famIsFamily(wave)
				;; call self with all family members
				foreach(sv famGetSweepValues(wave)
					wf = famValue(wave sv)
					icFamStatHop(wf work)
				)
			)
			( t
				error("icFamStat: cannot handle %L\n" wave)
			)
		); cond
	); let
)
;====================================================================

;====================================================================
;; icFamStat(type wave(s))
;; do statistic over a collection of waves/families
;====================================================================
procedure( icFamStat(type @rest waves)
	let( (work wf fm m v s k N scale)
		work = ncons(nil)
		
		;; run over all input waveforms
		foreach(wave waves
			icFamStatHop(wave work)
		)

		;; build result from collected data
		if( type=="non-central moments" then
			;; EX + non-central moments EX^2, EX^3, EX^4
			fm = famCreateFamily("statistics" 'string)
			wf = drCreateWaveform(work->xVec work->y1)
			famAddValue(fm "1st" wf)
			wf = drCreateWaveform(work->xVec work->y2)
			famAddValue(fm "2nd" wf)
			wf = drCreateWaveform(work->xVec work->y3)
			famAddValue(fm "3rd" wf)
			wf = drCreateWaveform(work->xVec work->y4)
			famAddValue(fm "4th" wf)

		else
			;; central/standardized moments
			;; sigma^2 = variance = E(X^2-u^2) = E(X^2) - E(X)^2
			;; skew     = E( ((X-u)/sigma)^3 ) = ( E(X^3) - 3*u*sigma^2 - u^3 ) / sigma^3
			;; kurtosis = E( ((X-u)/sigma)^4 ) = (E(X^4) - 4*u*sigma^3*skew - 6*u^2*sigma^2 - u^4) / sigma^4

			for(i 0 work->len-1

				; for unbiased variance: rescale from 1/N to 1/N-1 
				; better estimate for N>6, otherwise use biased variance
				N = work->N[i]
				scale = if(N<6 then 1.0 else N / float(N-1) )
				
				m = drGetElem(work->y1 i)
				v = drGetElem(work->y2 i)
				s = drGetElem(work->y3 i)
				k = drGetElem(work->y4 i)

				v = (v - m*m) * scale 
				if( v<=0 then
					v = 0
					s = 0
					k = 0
				else
					s = s - 3*m*v - m*m*m
					k = k - 4*m*s - 6*m*m*v - m*m*m*m
					
					unless( type=="central moments"
						v = sqrt(v)
						s = s / (v*v*v)
						k = k / (v*v*v*v)
					)
				)

				drSetElem(work->y2 i v)
				drSetElem(work->y3 i s)
				drSetElem(work->y4 i k)
			); for
			
			if( type=="central moments" then
				work->y2->name  = "variance"
				work->y3->name  = "u3"
				work->y4->name  = "u4"
			else
				;; standardized moments
				work->y2->name  = "stdev"
				work->y3->name  = "skew"
				work->y4->name  = "kurtosis"
			)
			
			fm = famCreateFamily("statistics" 'string)
			wf = drCreateWaveform(work->xVec work->y1)
			famAddValue(fm work->y1->name wf)
			wf = drCreateWaveform(work->xVec work->y2)
			famAddValue(fm work->y2->name wf)
			wf = drCreateWaveform(work->xVec work->y3)
			famAddValue(fm work->y3->name wf)
			wf = drCreateWaveform(work->xVec work->y4)
			famAddValue(fm work->y4->name wf)
		); if
		fm
	); let
)
;====================================================================


;====================================================================
;; GUI builder information
;====================================================================
ocnmRegGUIBuilder(
 '(nil
  function icFamStat
  name icFamStat
  description "Statistics over family graphs"
  category ("Custom Functions")
  analysis (nil
	general (nil
		args (type wave)
		signals (nil
			wave (nil
				prompt "Signal"
				tooltip "Signal"
			)
		)
		params(nil
			type (nil
				prompt "type"
				tooltip "statistics"
				guiRowHint 1
				default "standardized moments"
				type (
					"standardized moments"
					"central moments"
					"non-central moments"
				)
				required nil
			)
		)
		inputrange t
	)
  )
  outputs(result)
 )
)
;====================================================================

Hi Andrew:

The steps i'm following is like the follow,

1. plot the inlEndpoint(iout_sink), I'm running a DC simulation with vcode sweeping from 0 to 1023, and don't see it available in the Sweep Variable

2. plot swapSweep(inlEndpoint(iout_sink) "mcparamset" 500 )

3. plot stddev(swapSweep(inlEndpoint(iout_sink) "mcparamset" 500 )), The x-axis changed to temperature when I plot it, but the values are correct as it's taking the std deviation of code 500 for the 4 corners. I couldn't chage the x-axis back to vcode with the swapSweep as it seems like this variable doesn't exist.

Not sure if i'm doing it correctly, I don't understand the structure of the familywaveforms at the start, seems like there are a lot of parameters going on with corner simulation with monte carlo run.
For now henker's solution works ok with me so I'll use that solution for now, I'll go through the materials online to try and understand the various data structures for skill.

Thanks

Hi henker:

Thanks for the code, it works.
Just have a quick question on <list of waveforms or families> for icFamStat(<option> <list of waveforms or families>)

Currently I have 4 corner runs with 200 monte carlo runs, I only want to process the list of waveforms on C0_0 column shown in the red box.
However if I use icFamStat("standardized moments" inlEndpoint(iout_sink)) in the calculator, I feel like it's processing all the corners shown in the blue box.
Is there a way to know which kind of waveforms i'm processing? How do I force the waveform group in a single column like the red box? I'm not familiar with the difference of list of waveforms or families.
Of course a workaround is just to have 1 corner simulation for the monte carlo run to make sure it's getting the correct list of waveform and use the function but it's a bit tiresome.

Thanks for the reply. I wouldn't expect to see vcode in the sweep variables when you do the swapSweep, because vcode already is the x-axis and so there's no point changing it to what it already is. What's a bit strange though is that it's not showing up in the legend after that, so something seems awry.

Anyway, glad you've got a way forward with henker's solution (I didn't check henker's code in detail due to lack of time at the moment). Writing something specific in SKILL was always an option here - I was just trying to see if there was a way of doing this with built-in functionality. I'm always happy writing SKILL, but since I always say that you should check first if something exists already in the tools, I was trying to practice what I preach!

Andrew.

You can limit a multidimentional sweep  to a subset by using leafValue. Application is basically:
leafValue(<sweep bundle> <parameter name> <parameter value> ...)

You can find the actual documentation in the Virtuoso Visualization and Analysis XL manual.

However, the most easy way to get the proper naming is to plot your complete sweep, then rightclick a single wave in the waveform viewer (preferably one which you want to include in the subset) and choose "send to calculator". There you get the full expression to select this single wave out of the sweep. From this you can remove the unnecessary parameter restrictions to construct your desired subset.
Plotting this expression as is will also show you what waves the icFamStat will be working on when passed to it.

Additionally, you can also add the single waveforms one by one as list in the icFamStat call, so you can also calculate on parts of the actual sweep, but be aware that multiple entries of the same wave is not checked, these then just get higher weight in the calculation.